基于数值模拟的CO2地质封存项目井筒泄漏风险定量化评价方法

甘满光，雷宏武，张力为*，李小春，李 琦

(1.中国科学院武汉岩土力学研究所岩土力学与工程国家重点实验室，湖北武汉 430071；2.中国科学院大学，北京 100049)

CO2地质封存（GCS）是应对全球气候变暖的重要技术手段，其关键是将捕集的CO2安全长期地封存在深层地质构造（深部咸水层、无开采价值的煤层、枯竭油气藏等）中，防止CO2泄漏至浅层地下水层或逸入大气[1]。如果在封存过程中CO2发生泄漏，近地表动植物和地下水层均会受到不同程度危害[2–6]。联合国政府间气候变化专门委员会（IPCC）[7]将CO2泄漏的风险分为全球风险和局部风险。全球风险指封存的CO2泄漏到大气中，可能引发显著的气候变化。局部风险指CO2泄漏对泄漏点附近的地区环境、人体健康、动植物生长等产生影响。目前，关于CO2地质封存泄漏风险的相关研究大多是基于风险评价体系对CO2的泄漏风险和泄漏导致的环境危害进行评价[8–11]。

对CO2泄漏风险进行评价可分为定性评价和定量评价两大类。定性评价主要基于特征、事件和过程方法学来确定、分类和甄选影响CO2地质封存的各类因子，实现CO2泄漏风险的分级[12]。比较典型的有Oldenburg[13]开发的筛选和分级框架（SRF），该框架允许用户对CO2地质封存场地的地质特征属性进行任意权重和不确定性分配，对两个或多个站点进行风险等级相对评价和排序，并基于CO2泄漏引起的健康、安全和环境风险对潜在的CO2地质封存地点进行评估[14]。He等[15]基于贝叶斯网络（BN）开发了CO2封存的风险评估模型，介绍并分析了该模型在石炭系储层CO2注入过程的风险评价、CO2封存量预测和浅层地下水污染概率计算等方面的应用，并证实了模型应用于CO2封存风险评估的有效性。王重卿[16]通过分析CO2地质储存安全风险因素，从地下和地面两方面建立层次分析结构的安全风险评价指标体系并确定各项指标权重，对CO2地质封存安全风险进行了综合评价。

定量评价指利用科学、系统的风险评估工具来预测CO2的泄漏量，进而采取有效的预防或补救措施，使CO2泄漏风险最小化，是评估CO2泄漏风险的重要手段[17–18]。胡丽莎等[19]将CO2泄漏风险源分为点源风险源、线源风险源及面源风险源，并给出了不同泄漏情景下CO2地质封存项目环境风险源强的计算方法。Oxand公司开发的P&R软件[20]可评价水泥降解和套管腐蚀导致的CO2沿井筒泄漏的风险，该软件应用于德国的Ketzin项目[21]，绘制了CO2在井筒内的路径分布和CO2泄漏风险分布图，对CO2的泄漏风险水平进行了量化和排序。此外，Tanaka等[22]开发的CO2地质封存环境风险评估系统可用于评估CO2泄漏量及CO2泄漏对注入点附近地面、海洋环境和周围空气的影响。

由于CO2地质封存系统的复杂性，CO2在储层中的分布和CO2泄漏通道的位置、数量和井筒渗透率等均存在不确定性，直接应用于定量评价CO2泄漏存在困难[23]。现有量化评价方法主要存在的问题表现在[24–25]：1）井筒和储层的控制方程不同，需要在统一框架内求解；2）井筒和储层在时间和空间尺度上存在显著差异，需要有效、可行的偏微分方程组进行求解；3）井筒与储层之间相互作用导致的其他复杂性。本文介绍了自主研发的可量化评估CO2泄漏风险的WellRisk软件，来定量预测中国首个CO2咸水层封存示范工程——神华鄂尔多斯CO2咸水层封存场地CO2沿井筒（注入井和监测井）的泄漏量，并将模拟结果与NRAP–IAM–CS的模拟结果进行对比，评估了CO2沿井筒的泄漏风险。井筒和储层的流体由完全不同的物理场控制，特别是井筒中的两相非达西流动和储层中的两相多孔介质达西流动有本质区别。为解决上述问题，WellRisk软件采用了一种模拟井筒–储层系统中非等温、两相和多组分流动的方法，将井筒和储层中的流动和运输作为一个集成系统求解，基于混合隐式方程，通过漂移流模型求解井筒中的瞬态动量方程，使得井筒中的流动与储层中的达西流无缝耦合。由于CO2泄漏过程的不确定性因素较多，本文的量化分析过程选取了对CO2泄漏影响较大的泄漏系数作为分析参数，得到不同泄漏系数条件下CO2沿井筒的泄漏量，以期评估神华鄂尔多斯CO2咸水层封存项目长期封存的安全性。

1 基于漂移流模型的CO2泄漏计算软件开发

CO2地质封存场地的井筒系统是CO2泄漏的主要通道，如图1所示[26]。其中，包括围岩–水泥环交界面、水泥环–套管交界面以及井筒水泥内部的孔裂隙通道。为保证CO2地质封存项目正常运行，需保证注入的CO2不发生泄漏。当封存的超临界CO2发生泄漏时，随着地层温度和压力的降低，CO2将在较浅的地层经历超临界状态转变为亚临界状态的过程，并逐渐气化逸出地面[27]。CO2相态的转变将引起其物性发生剧烈变化，从而对CO2地质封存场地的CO2泄漏过程产生影响[28–29]。Cai等[30–31]以神华鄂尔多斯CO2咸水层封存场地为研究对象，对注入过程中井筒–储层系统内的流体行为进行了分析。结果表明，在间歇注入过程中的停注期，井筒温度在围岩热传导作用下上升，而压力在储层调节下降低，导致300m以上深度井筒内的液相CO2发生气化，井筒–储层系统中的液相或超临界CO2发生气化后将促进CO2的泄漏。

图1 CO2地质封存场地井筒系统潜在CO2泄漏通道[26]Fig.1 Potential CO2 leakage channel of wellbore system at CO2 geological storage site[26]

CO2地质封存过程中，与井筒泄漏相关的CO2在井筒环空中的流动问题研究较少，主要是由于井筒环空具有很高的孔隙率和渗透率，流体在其中的流动属于非达西流范畴。传统的井流模型没有考虑相间的速度差异，认为流体在井中以相同的速度流动，这样会导致流体分数和压降预测结果不准确。多孔介质中的流动问题可以用传统达西定律求解，对于井中的流动则需要找到合适的描述方程，如漂移流模型[32–33]。漂移流模型的基本思想是将气液两相看作由气液性质差异导致相滑移的一个有效流体相，用混合速度表征气液各相的速度[24,34]。该模型考虑了流体相间的速度差异，可以模拟流体在不同速度下的流动情况以及瞬时流问题，是模拟CO2沿井筒泄漏的理想方法。考虑到井筒中多相流体的快速流动过程，对传统的漂移流模型进行了改善，用瞬时动量平衡方程代替稳态压力损失方程来求解混合速度问题，以提高模拟结果的精确性[34]。在求解方法上，漂移流模型涉及的质量、动量和能量守恒方程均可用数值方法求解，井中的热传导过程则用半解析解法或数值解法求解。在求解过程中，组分质量和能量平衡方程的空间离散用传统的积分差分法，时间的离散用反向1阶全隐式有限差分法。漂移流模型用牛顿迭代法对方程进行计算求解，在计算过程中，由于井流速度较快以及井流对时间的敏感性，需要设置较小的时间步长[34]。

基于上述漂移流模型，中国科学院武汉岩土力学研究所开发了国内首个CO2地质封存条件下封存CO2沿井筒泄漏量的WellRisk计算软件，通过导入储层CO2饱和度、压力、温度分布的模拟数据，计算与CO2储层交汇井筒中的CO2泄漏量，从而为CO2地质封存工程的长期封存可靠性和环境风险评价提供数据支撑。WellRisk软件旨在满足国内CO2地质封存示范项目对于CO2泄漏风险定量评估的现实需求，实现模拟结果与CO2泄漏风险评价软件NRAP–IAM–CS模拟结果的相互印证，并突破NRAP–ICM–CS软件具有的无法考虑CO2沿井筒泄漏过程中相态变化的局限。软件的基本架构如图2所示。

图2 自主开发的WellRisk软件的基本架构Fig.2 Basic structure of the self-developed WellRisk software

2 应用实例

2.1 神华鄂尔多斯CO2咸水层封存项目

神华CCS（CO2capture and storage）项目是中国CO2捕集与地质封存的首次规模化探索，是中国首个全流程咸水层CCS示范工程[35]。CO2主要来源于煤液化和煤制烯烃和甲醇等燃料的过程，项目从2011年开始注入至2015年封井，注入CO2约30×104t，目前处于注入后监测阶段。神华CO2地质封存场地注入井和监测井的相对分布以及地层结构如图3所示[36–38]。项目位于内蒙古自治区鄂尔多斯市东南部，CO2注入场地距离捕集场地约17 km，场地内布置了一口注入井（INJW）和两口监测井（MW1和MW2）。

图3中，监测井 1 位于注入井正西70 m，负责监测不同层位的孔隙压力和储层温度；监测井2位于注入井正北30 m，主要用于盖层以上的含水层中取水样以及进行时移垂直地震剖面（VSP）测井[37]。场地内发育的地层从上往下依次有三叠系纸坊组、和尚沟组、刘家沟组，二叠系石千峰组、石盒子组、山西组，奥陶系马家沟组。通过对地层的孔隙率、渗透率及矿物成分分析，将上至刘家沟下至马家沟在内的1 576～2453m地层作为注入储层[39]。

2.2 神华CCS场地CO2沿井筒泄漏量计算

WellRisk软件采用的漂移通量法将CO2–咸水两相流看作一个整体，整体具有一个混合速度，但其内部由于CO2和咸水的性质差异，存在相滑移，且CO2和咸水两相各自存在一个流动速度（是混合速度的函数）。当CO2从井筒中发生泄漏时，其相态会由超临界态转变为气态，这个过程会导致CO2的密度和黏度等性质发生变化，进而导致CO2的流动速度发生变化，而这些变化过程在漂移通量法的控制方程中均可以体现。漂移流模型假设气体速度uG通过经验本构关系与两相流混合流速j和气体的漂移速度ud相关联，所有速度变量在截面上取平均值[24-25]，具体见式（1）：

式中，C0为反应流动状态系数（剖面参数），用于考虑局部含气饱和度和速度剖面对管道截面的影响，其中紊流时取1.2，层流时取2.0。

两相流混合流速j和咸水流速uL可分别写为式（2）和（3）：

式（2）～（3）中，SG为CO2的气相饱和度，井筒内两相流动量方程可简化为混合速度um和漂移速度ud的单一方程，具体可见式（4）[24–25]：

式中： Γ为横截面周长；f为壁面剪应力，是壁面粗糙度和流动状态的函数；A为随z（沿井筒的空间坐标）改变的井筒横截面面积； γ为表征气液两相滑移程度的系数， γ的计算过程见式（5）：

其中：

式中， ρ为密度，S为饱和度，u为流速，P为压力，Se为CO2发生泄漏的有效截面积， θ为井筒截面与垂直方向的夹角，下标G、L为流体相位（气体和液体），m为气液混合物，t为时间，其余变量的含义同上。

计算CO2沿井筒泄漏量的关键是确定CO2发生泄漏的有效截面积（Se）。对于开放井筒，Se即为井筒的总截面积。然而，在实际工程中，由于封隔器、井筒水泥的存在及严格规范的施工流程保证了固井和封堵的质量，因此，几乎不存在CO2沿井筒总截面泄漏的可能。一般情况下，CO2只能通过固井水泥–套管交界面、固井水泥–围岩交界面、油管–套管环空等途径发生泄漏，对应的Se值远小于井筒的总截面积。总体而言，Se值与井筒的固井质量密切相关，固井质量较好的井筒，其对应的Se值较低。在现场操作中，一般通过井筒完整性测井确定固井质量，但由井筒完整性测井曲线获得有效截面积仍存在困难，典型的井筒完整性测井曲线只能定性回答井筒固井质量是否良好，难以直接推导出Se的数值。因此，本文给出了一个简便易行的解决方案，即采用式（10）计算Se：

式中：ST为井筒总截面积；k为泄漏系数，即不同固井质量的井筒截面有效泄漏面积与截面面积的比值，根据固井质量确定。本文分别测试了泄漏系数为10–5、10–6、10–8的情况，基于模拟结果，给出了上述取值对应的井筒固井质量。

如图4所示，利用WellRisk软件计算CO2沿井筒泄漏量的一般流程为：首先确定拟进行研究的井筒，为井筒选择对应的井筒泄漏系数，计算井筒的有效环空截面积Se；然后，导入井筒所在位置的储层CO2饱和度、压力、温度分布；最后，在软件的用户界面输入井筒深度、粗糙系数、换热系数等参数，运行软件，即可获得CO2沿井筒的泄漏量。计算过程基于校准的储层模型，采用ECLIPSE 300软件模拟获得不同时间点的30组储层数据（CO2饱和度、温度、压力等）作为输入数据，每1组储层数据对应相隔一定时间内（模拟总时间为10 0 0 a）注入井与监测井泄漏量的模拟结果，即基于ECLIPSE 300的计算结果，储层的数据随时间（CO2注入期间和注入后阶段）发生变化，储层数据在WellRisk软件中是随时间更新的。

图4 WellRisk软件几何计算模型（输入与输出参数）Fig. 4 Geometric calculation models (input and output parameters) for WellRisk

根据神华CCS现场场地参数，本次模拟采用的相关参数见表1。

此外，为了对WellRisk的模拟结果进行验证，采用了美国能源部开发的NRAP–IAM–CS软件开展相关的数值模拟，模拟的具体流程如下：

在NRAP–IAM–CS软件中，CO2储层参数（储层深度、温度、压力、CO2饱和度、溶解CO2浓度以及渗透率）由相关油藏/地球化学数值模拟软件，如FEHM[41]、TOUGH2[42]、ECLIPSE[43]、CMG-GEM[44]等。根据CO2地质封存场地的实际地层条件模拟得到。根据神华CCS地质封存场地的实际地层条件，首先使用模拟软件ECLIPSE 300[45–46]模拟CO2的注入过程，将注入井深度1631m处压力的模拟结果与场地监测数据进行对比以验证模型的可靠性[46]。然后，通过改变地层岩石/流体参数，得到不同岩层和流体参数条件下储层参数的空间分布。为解决有限储层数据的不确定性，基于校准的储层模型[47]和相关CO2地质封存现场的实际场地调研[48]，构建了30组储层模拟数据，以涵盖砂岩渗透率、相对渗透率、毛细管压力和盖层渗透率等一系列参数的不确定性，具体数值参见文献[49]。

表1 W ellRisk软件的输入模拟参数Tab.1 Input simulation parameters for WellRisk software

浅层地下水的物理、水文地质和地球化学参数以及储层和含水层的基本物理参数参见文献[49]和表2。

表2 井筒泄漏简化模型的输入参数Tab.2 Input parameters for the simplified model of wellbore leakage

参照神华CCS项目的场地条件，算例中整个系统包括3个井筒（1个注入井和2个监测井），井筒的位置根据神华CCS项目实际场地的井筒相对位置而定[38–39]。根据Carey[50]的调研，渗透率高于10–12m2的井筒被认为是高泄漏风险井，因此，将神华CCS现场注入井和监测井固井质量很差的井筒渗透率设定为10–11m2、固井质量较差的井筒设定为10–12m2和固井质量良好的井筒设定为10–14m2。根据神华CCS项目的实际注入情况，将软件运行时间设定为1000a（3a注入+997 a注入后阶段），以计算CO2沿井筒的长期泄漏量。

计算CO2沿井筒的泄漏量时，将储层内CO2的泄漏视为多孔介质的流动，根据不同几何形状泄漏通道和渗透率之间的关系，对CO2的泄漏量进行计算[51]，如式（11）：

式中，Q为体积流量，keff为井筒渗透率，A为井筒横截面积，ΨL为泄漏通道源处的势，ΨT为泄漏通道顶部处的势，L为泄漏通道的长度。与W ellRisk软件不同，NRAP–IAM–CS计算公式中不含有CO2密度相，因此，NRAP–IAM–CS无法考虑CO2发生相变时导致的密度显著变化。

3 结果与讨论

图5为利用WellRisk软件运行的CO2沿泄漏系数10–6井筒的累计泄漏量变化。由图5可知，CO2沿井筒的泄漏量随时间基本呈现线性变化。图6为W ell-Risk软件预测的CO2沿泄漏系数10–6和10–8井筒（注入井和监测井）的泄漏量，CO2沿泄漏系数为10–6井筒的泄漏量明显大于泄漏系数10–8井筒的泄漏量，表明泄漏系数k越低，CO2的泄漏量越少。值得注意的是，当泄漏系数为10–5时，井筒的泄漏量预测结果出现异常（泄漏量大于注入量）。出现这种情况的原因是，WellRisk无法实现泄漏对储层饱和度的逆反馈，亦即大量CO2泄漏不会影响储层中的CO2饱和度输入数据，导致CO2泄漏量很大时，总的泄漏量甚至可能大于注入总量。在此种情况下，具体的泄漏数值不具实际意义，但此种情况表明CO2存在严重泄漏。

图5 CO2沿泄漏系数为10–6井筒的泄漏量变化（1000a）Fig.5 Changes in leakage of CO2 along a wellbore with an annular coefficient of 10–6(1000a)

图6 CO 2分别沿泄漏系数10–6和10–8的井筒泄漏量对比（1000a）Fig.6 Comparison of total CO2 leakage along the wellbore with leaakage coefficients of 10–6 and 10–8(1000a)

CO2泄漏至大气的量随井筒渗透率的降低大幅减少，将CO2泄漏到大气的量除以神华CCS场地总的CO2注入量（0.3×106t），可得到CO2的泄漏百分比，该指标是评价CO2地质封存场地CO2封存安全性的关键指标。根据IPCC的报告[7]，对CO2地质封存工程的现场监测结果表明，在安全的地质储层中CO2泄漏占CO2注入总量的比例应在1000a内低于1%。因此，选择1000a内注入CO2总量的1%作为评估CO2泄漏风险的阈值。CO2沿泄漏系数10–6和10–8井筒（注入井和监测井）的总泄漏量分别为7.20×102t和5.49t，经计算可得，1000a后CO2经泄漏系数为10–6、10–8的井筒泄漏至大气的泄漏百分比分别为0.24%、0.002%，井筒泄漏百分比低于IPCC建议的风险阈值。因此，泄漏系数为10–6、10–8的井筒属于低泄漏风险井筒，对应较好的固井质量。而泄漏系数为10–5井筒的泄漏量巨大，因此，泄漏系数为10–5的井筒属于高泄漏风险井筒，对应较差的固井质量。

图7为利用NRAP–IAM–CS软件计算的CO2沿不同渗透率井筒向大气泄漏量随时间的变化以及按平均值计算的1000a后总的CO2泄漏量（假定注入井和监测井的渗透率分别为10–11、10–12和10–14m2）。图7中，10%、50%、90%和平均分别表示概率为10%、50%、90%和平均的CO2泄漏量。由图7可以看出，CO2沿渗透率为10–11、10–12和10–14m2的井筒泄漏至大气中的CO2量随时间线性增大，表明注入结束后CO2沿井筒的泄漏速率基本不变。其中，CO2沿渗透率为10–11m2井筒的泄漏量与为10–12m2井筒的泄漏量为同一数量级，比渗透率为10–14m2的井筒的泄漏量高一个数量级。

经计算可得，CO2经渗透率10–11、10–12、10–14m2的井筒泄漏至大气的量分别为5.86×103、1.46×103和9.00×102t，算得泄漏百分比分别为1.95%、0.49%、0.30%。因此，渗透率为10–11m2的井筒泄漏百分比高于IPCC设定的安全标准，是高泄漏风险井筒，而渗透率为10–12、10–14m2的井筒1000a后的泄漏百分比小于1%，说明这两个渗透率的井筒是安全的。此外，Carey[50]指出，渗透率大于10–12m2的井筒通常具有贯通环空带、水泥塞封固质量较差等明显影响井筒封固性能的缺陷，因此是有泄漏风险的井筒。模拟结果较好地验证了Carey的结论。

值得注意的是，图7所示的CO2泄漏情形与现实存在一定的差异，这是由于NRAP–IAM–CS软件中无法考虑CO2沿井筒泄漏导致的储层CO2饱和度降低的行为，即储层的CO2饱和度变化忽略了CO2沿井筒的泄漏情形。而在实际的封存场地中，储层的CO2饱和度随着CO2沿井筒的泄漏持续降低。因此，图7中的结果代表了CO2泄漏的最坏情况。神华CO2地质封存场地的CO2注入井和监测井井龄均在10 a左右，固井质量良好，从测井数据看无明显缺陷[52]。因此，可以预测神华项目CO2通过注入井和监测井泄漏的百分比低于IPCC的安全阈值，该工程是安全的。从现场的3次VSP监测数据也可以证实模拟结果[53–54]，现场监测结果表明，注入CO2的运移停留在注入井周围450m以内，没有CO2从井筒垂直泄漏到地表，这进一步证明了神华CO2地质封存场地进行CO2封存的安 全性。

图7 CO 2沿不同渗透率井筒的泄漏量随时间的变化以及按平均值计算的1 0 0 0 a后CO 2沿不同渗透率井筒的总泄漏量对比Fig.7 Comparison of CO2 leakage along wellbores with different permeability over time and the total leakage of CO2 along wellbores with different permeability after1000 years(calculated based on average)

图8 神华CO 2地质封存场地注入点附近地下水pH值的变化[53]Fig.8 Change of pH value of groundwater near the injection point of Shenhua CO2 geological storage site[53]

此外，Zhao等[53]对神华CO2地质封存场地注入点附近的地下水的pH值进行了监测，结果如图8所示。由图8可见，从2011年1月至2014年8月，地下水的pH值未发生降低，由于CO2泄漏会与地下水形成碳酸，导致地下水的pH值减小，因此该结果表明无CO2从储层发生泄漏。李琦等[55–56]利用U型取样管在2015年5月、8月和10月对地下16m和2 0 m的地下水进行了取样分析，结果发现，20m深地下水的pH值和Cl–、Ca2+、Mg2+等离子浓度的水平低于16m的情况。中国地质调查局对该地区的监测数据显示[53]，注入点附近地下水的pH值和矿化指标随时间变化不明显，因此，监测数据的差异在正常范围内，表明CO2没有泄漏至地下水。

通过对比WellRisk和NRAP–IAM–CS软件的计算结果可以看出：在WellRisk软件中，当井筒泄漏系数设定为10–5时，模拟结果远远高于NRAP–IAM–CS的模拟结果，说明WellRisk软件在高泄漏井筒模拟方面存在误差；当井筒泄漏系数设定为10–6时，模拟得到的CO2沿井筒的泄漏量与NRAP–IAM–CS中CO2沿渗透率为10–14m2井筒泄漏的结果相对接近；当井筒泄漏系数设定为10–8时，WellRisk的模拟结果比NRAP–IAM–CS的结果小两个数量级，此条件下泄漏量绝对值本身很低，即使按较高的NRAP–IAM–CS的模拟结果，1 000 a后的CO2泄漏量仅占总注入量的0.3%，故两个模拟结果均表明此种情况下基本无CO2泄漏，CO2泄漏风险极低。

从模拟结果也可以看出，WellRisk软件的漂移流模型用混合速度表征气液各相的速度，利用有效流体相代表气液两相，简化了模型，但同时也增加了不确定性。在WellRisk软件中，CO2的泄漏不影响储层的CO2饱和度输入数据，选定时间内由ECLIPSE软件计算的储层输入数据无变化。实际上，当CO2发生泄漏后，储层的CO2饱和度会降低，导致CO2的泄漏量会进一步降低。因此，WellRisk软件简化了储层CO2的泄漏过程，得到的泄漏量偏大于实际CO2地质封存场地CO2的泄漏量。综上所述，WellRisk软件模拟的是CO2地质封存场地最不利的CO2沿井筒泄漏情形，这对决策者分析泄漏风险的最不利情况是很有帮助的。此外，根据模拟结果可知，静态的井筒渗透率是可预测的，在CO2地质封存的高温高压条件下，关键是确定CO2水岩反应对井筒渗透率的影响，进而确定CO2的泄漏风险，由CO2水岩反应导致的井筒渗透率动态变化的相关模块目前正在开发。

4 结 论

本文采用自主研发的可定量预测CO2地质封存场地的CO2沿井筒泄漏量的WellRisk模拟软件预测了中国首个CO2咸水层封存项目，即神华CO2地质封存项目的CO2沿井筒的泄漏量，定义了WellRisk软件运行所需的井筒泄漏系数，并与已有的NRAP–IAM–CS软件对CO2沿不同渗透率井筒泄漏量的结果进行了对比，评估了CO2沿井筒泄漏的风险，得到以下结论：

1）CO2沿井筒泄漏至大气的量随着井筒泄漏系数的降低而减少，泄漏系数为10–5井筒的CO2泄漏量巨大，远大于风险阈值，因此是高风险泄漏井，而泄漏系数为10–6、10–8井筒的CO2泄漏量占注入总量百分比小于风险阈值，因此是低风险泄漏井。基于上述结果，本文建议当井筒固井质量良好时，泄漏系数取10–6或更低，而当井筒固井质量较差时，泄漏系数取10–5或更高。

2）CO2沿井筒泄漏至大气的量随着井筒渗透率的降低而减少，渗透率为10–11m2井筒的CO2泄漏量占注入总量百分比大于风险阈值，因此是高风险泄漏井，而渗透率为10–12、10–14m2井筒的CO2泄漏量占注入总量百分比小于安全阈值，因此是低风险泄漏井。此结果与之前Carey[50]报道的渗透率高于10–12m2的井筒与高泄漏风险井筒的调研结果相吻合。

3）神华场地的井筒使用年限较短，固井效果良好，对应的井筒泄漏系数为10–6或更低，渗透率应为10–14m2或更低。根据本文获得的模拟结果可以预测，神华场地的CO2泄漏量很少，泄漏风险低，CO2可长期安全封存在地下，现场的CO2泄漏监测结果也证实了神华场地的封存安全性。

综上所述，本文根据神华CO2地质封存场地储层CO2饱和度分布、压力分布和温度分布数值模拟结果，利用WellRisk软件计算了CO2沿不同泄漏系数井筒的泄漏量，并与NRAP–IAM–CS软件的结果进行了对比，验证了自主研发的WellRisk软件的可靠性。泄漏系数和渗透率是两种不同计算方法分别采用的参数，泄漏系数的定义在本文中首次提出，而NRAP–IAM–CS软件采用的大渗透系数法的控制方程是达西定律的形式，渗透率是该方程中的参数，用于描述流体在井筒内流动的难易程度，与井筒截面有效泄漏面积无关。因此，泄漏系数和渗透率无直接的关系。本文首次提出并定义了井筒泄漏系数，结合NRAP–IAM–CS软件的泄漏风险定量计算结果，为不同的井筒泄漏风险提供了相对应的井筒泄漏系数取值。本文实现了考虑井筒泄漏系数和井筒渗透率的CO2沿封存场地井筒泄漏量的定量计算，预测了CO2沿不同泄漏系数和渗透率井筒的泄漏风险，对CO2地质封存场地的长期安全性进行了评估与验证。