原地生成宇宙成因核素暴露和埋藏年齡計算方法的統一*

黃費新 程楊 李巖 李廣偉 董國成 梁霞

1. 中國冶金地質總局礦產資源研究院，北京 101300 2. 南京大學地球科學與工程學院，南京 210023 3. 中國科學院地球環境研究所，西安 710061 4. 西安加速質譜中心，西安 710061 5. 中國地質科學院地質力學研究所，北京 100081

原地生成宇宙成因測年技術中，暴露年齡、埋藏年齡及侵蝕速率是其中三個重要研究方向(Dunai, 2010; Gosse and Phillips, 2001; Lal, 1991; Granger and Muzikar, 2001; 孔屏, 2002, 2012; Blardetal., 2019; 黃費新等, 2019)。其中暴露年齡的研究得到了最為廣泛的應用，國內的相關研究也很多(黃費新等, 2004, 2005; 李廣偉等, 2009; 董國成等, 2014; Kongetal., 2010, 2011; 張金玉等, 2018)。埋藏年齡的研究方法也逐漸在國內得到推廣(Kongetal., 2012; Liuetal., 2013; 劉彧等, 2013; 鄭勇等, 2014; 韓非等, 2016; 陳清敏等, 2018)。侵蝕速率方面，受方法學本身的局限，大多僅限于假設侵蝕速率恒定時計算得到最大侵蝕速率(Lal, 1991; Kongetal., 2007)或侵蝕速率有限次數變化的情況(Carracedoetal., 2019; Granger and Riebe, 2014)，對侵蝕速率連續變化時的研究探索非常少(Knudsenetal., 2019)。其中，計算埋藏年齡的關鍵，應是埋藏理論模式的建立，從而確定埋藏前核素的濃度，以及埋藏過程中核素有無新生量。目前，埋藏理論模式中快速深埋藏模式比較受到重視，這種忽略埋藏速率大小，把埋藏過程看成極短時間完成的最簡單化的埋藏模式，除了洞穴沉積物比較接近實際情況外，其它情況可能相距實際較遠從而帶來很大的研究誤差(Knudsen and Egholm, 2018)。自然界較普遍的情況，應是具有埋藏前濃度的樣品，埋藏深度在埋藏速率變化的情況下發生變化(若埋藏深度不變視為埋藏速率為0，埋藏深度減小視為埋藏速率為負值)。根據埋藏前樣品是否具有原始濃度，以及埋藏速率的不同變化情況，可以建立不同的埋藏理論模式，并推導出不同埋藏理論模式下核素濃度的計算等式，進而計算埋藏年齡。暴露和埋藏在以往的研究中是當成兩個彼此相對獨立的研究領域，并且彼此間往往成為研究中的干擾因素。雖然已經有研究者開始嘗試推導針對侵蝕和埋藏過程利用拉格郎日方法(以樣品為研究關注點)或歐拉方法(以地表位置為關注點)推導核素濃度的數學計算方法(Knudsenetal., 2019)，但尚沒有明確提出將侵蝕速率和埋藏速率視為同一參數。實際上，如果將侵蝕和埋藏對樣品距離地表深度的改變效果關聯起來考慮，它們其實是一種反向過程：暴露或埋藏過程中，侵蝕使樣品距離地表越來越近，埋藏使樣品距離地表越來越遠。

原地生成宇宙成因放射性核素的濃度和暴露年齡、侵蝕速率之間的通用計算等式為(Lal, 1991)：

(1)

(1)式為樣品有原始濃度時的核素濃度與暴露時間的計算關系式，等式右邊第一項為核素新生濃度，第二項為暴露年齡計時前已有的濃度衰變后的繼承濃度。式中N是目前樣品中的放射性核素濃度(atoms/g)，N0是暴露年齡計時前已有的濃度。P是核素的地表生成速率(atoms/(g·yr))，數值隨地磁緯度及高程的變化而變化。λ是核素的衰變常數(穩定核素可視為λ=0的放射性核素)。ρ是巖石密度(g/cm3)。ε是地表侵蝕速率(cm/yr)。Λ是宇宙成因核素的吸收深度(g/cm2)。T是樣品的暴露年齡(yr)。

如果暴露計時前樣品無原始濃度即繼承濃度為0，得：

(2)

(3)

如果侵蝕速率很小可忽略，(2)式改寫為得到廣泛運用的最小暴露年齡的計算等式：

(4)

等式(1)中繼承濃度這一項N0e-λT也能代表樣品中的原始放射性核素濃度隨暴露時間增加而衰變的過程。埋藏年齡的計算方法，與這一衰變過程的濃度變化計算原理密切相關。最簡單的情況，即在快速深埋藏模式下，設樣品在暴露期間的最后濃度是N0，其后迅速埋藏到足夠深度(通常認為大于20m)，不再產生新的核素，只有衰變過程，則樣品中的核素濃度和埋藏年齡(T2)的計算關系如下(注意此時等式中代表原始濃度的N0較(1)式已經有新的地質含義，即代表的已經是埋藏前已有濃度，而不是暴露計時前已有的濃度)：

N=N0e-λT2

(5)

(5)式其實就是利用了放射性核素濃度衰變計算等式，是計算快速深埋藏年齡的實質性核心原理，求解過程則只不過是計算技巧和參數代入方面的問題：利用兩種放射性核素(比如常用的10Be和26Al)，將其埋藏前的原始濃度N0替換為(3)式即假設埋藏前為穩態侵蝕(增加未知數埋藏前侵蝕速率ε)或(4)式即假設0侵蝕速率(增加未知數埋藏前暴露時間T1)的濃度計算關系式，就可聯立方程解出兩個未知數T2和ε(或T1)(Granger and Muzikar, 2001; 孔屏等, 2012; Granger, 2014)。有時也采取更簡單的策略，即假設埋藏前的核素濃度比值就等于核素生成速率比值，埋藏年齡就是由等于核素生成速率比值的原始濃度比值因衰變而逐漸減小，最后達到實測核素濃度比值所需要的時間。當然只有在埋藏前暴露時間較短或者侵蝕速率很高時，樣品的埋藏前核素濃度比值才接近其生成速率比值(Lal, 1991; Matmonetal., 2014)。由于放射性核素的濃度衰變計算等式即(5)式在溫度、壓力等物理條件變化下是不變的，所以，求解不同埋藏理論模式下的埋藏年齡，關鍵是要確定出每種核素在該種埋藏模式下的核素原始濃度、埋藏后有無新生濃度，以及有新生濃度時的濃度生成規律，然后求解埋藏過程涉及到多少個獨立未知數，就要測試多少種放射性核素(不同穩定核素只能建立一個同質的濃度計算等式，僅能視為同一種放射性核素)，聯立多少個濃度方程來求解。求解方程的過程其實不應該是重點，如果不能得到代數解析解，至少是可以采用賦值逼近法計算的。

計算最小暴露年齡的(4)式和描述快速埋藏過程的(5)式，其實都是可統一于通用等式(1)中的。其它較為復雜的埋藏過程，是否也能統一于等式(1)中呢？為了探討這個問題，并建立不同埋藏理論模式下的濃度計算等式，首先需要探討埋藏和侵蝕在核素濃度計算中所起作用的區別和聯系。

1 埋藏速率可視為負侵蝕速率的數學證明

宇宙成因核素的生成速率可劃分為三個主要部分，即散裂反應、快μ介子和負μ介子對應的生成速率部分(Dunai, 2010; Heisingeretal., 2002a, b)。散裂反應對應的生成速率部分，隨樣品距地表深度x(cm)的增加呈指數關系快速下降(即P樣=P地表e-ρx/Λ散)，快μ介子和負μ介子對應的生成速率隨深度增加基本上也是呈指數關系下降的，但下降速率要慢，即散裂反應、快μ介子和負μ介子對應的核素吸收深度Λ數值不同。由于在地表位置，快μ介子和負μ介子對地表核素生成速率貢獻比例不大，在計算地表樣品的暴露年齡時往往不區分散裂反應、快μ介子和負μ介子對應的生成速率部分，而綜合成地表生成速率。在計算埋藏年齡時，由于快μ介子和負μ介子對應的生成速率隨埋藏深度的增加而下降的速率相對較慢，隨深度增加在核素生成濃度的占比中越來越高。埋藏至一定深度(一般認為20m)后，散裂反應、快μ介子和負μ介子對應的生成速率都將下降到0。除了快速深埋藏模式，為減小誤差，淺埋藏、慢速深埋藏模式計算核素濃度與埋藏時間關系時，應該區分散裂反應、快μ介子和負μ介子對應的生成速率部分對核素總濃度的影響。

針對有侵蝕的暴露過程，先僅考慮散裂反應對應的生成速率部分，xcm深度處無原始濃度的樣品濃度變化的微分方程是dN/dt=-λN+P樣=-λN+Pe-ρx/Λ(此處P專指樣品對應的地表散裂反應核素生成速率，不是樣品所在深度xcm處的散裂反應核素生成速率，Λ專指散裂反應對應的吸收深度)。設總暴露時間為T，侵蝕速率ε恒定，在0-T時間內，隨著侵蝕的進行樣品從深部逐漸接近地表，t時刻樣品距地表深度x=ε(T-t)，濃度變化的微分方程將是dN/dt=-λN+Pe-ρε(T-t)/Λ，解此微分方程得到等式(2)(黃費新等, 2019)。注意所研究的樣品因侵蝕逐漸變淺，T時刻即目前已經在地表，樣品的散裂反應核素生成速率P樣已經等于地表散裂反應核素生成速率P。

埋藏過程中，設無原始濃度的樣品初始時刻在地表，在埋藏過程持續的總時間T內以恒定埋藏速率β(burial第一個字母)逐漸埋藏于地表之下(注意埋藏過程中也有侵蝕存在，但可視為抵扣了最上層的埋藏量，抵扣了埋藏過程中的侵蝕量后的凈埋藏速率才是β，因此侵蝕量在推導埋藏過程中濃度的計算等式時是已經抵扣過的)。那么埋藏過程中的t時刻樣品的深度x=βt。此時的濃度變化微分方程為：

(6)

此微分方程解法完全類似dN/dt=-λN+Pe-ρε(T-t)/Λ且更簡單(黃費新等，2019)，解得：

(7)

(7)式中P仍為地表核素生成速率，但當前時刻樣品已經在地表下βTcm處(即采樣深度xcm)，故其生成速率為P樣=Pe-ρβT/Λ，代入(7)式得：

(8)

比較(2)式和(8)式，形式幾乎一樣，差別處是(2)式的P在(8)式中為P樣，λ+ρε/Λ換為λ-ρβ/Λ，相差一個正負號。(2)式中由于地表核素生成速率即樣品的生成速率，因此(2)式和(8)式可以統一，意味著以樣品當前所在深度的生成速率為參照，勻速埋藏速率可視為負勻速侵蝕速率。樣品的濃度計算等式統一為：

(9)

那么，(9)式中參數ε既可代表侵蝕速率，又可代表埋藏速率(視為侵蝕速率取負值，即β=-ε)，侵蝕速率和埋藏速率可以統一由同一參數指代。P樣是樣品在目前位置處(核心要素是樣品距地表深度和所在地質體密度)的核素生成速率。

由于核素衰變快慢與樣品所處深度無關，同理，等式(1)中將P替換為P樣，再視ε=-β，就不僅是勻速侵蝕時暴露年齡的計算等式，也是勻速埋藏時埋藏年齡的計算等式。即：

(10)

快μ介子和負μ介子對應的生成速率部分，與散裂反應對應的生成速率部分隨樣品距地表深度下降呈指數函數變小的規律是相同的，同理可知埋藏速率可視為負侵蝕速率的原理對快μ介子和負μ介子對應的生成速率部分也是成立的。

非勻速侵蝕和非勻速埋藏下，可以利用(8)式和(9)式分別對t求導，得到濃度瞬時變化的計算式，前者N′=P樣e-(λ-ρβ/Λ)t，后者N′=P樣e-(λ+ρε/Λ)t，埋藏速率仍然可視為負侵蝕速率，由于每一瞬間埋藏速率都相應可視為負侵蝕速率，最終求積分結果埋藏速率可視為負侵蝕速率的規律仍然成立。由于侵蝕速率或埋藏速率連續變化的情況下，求濃度的微積分過程中涉及到形如e-t2函數的求積分，其結果不是初等函數，所以侵蝕速率或埋藏速率連續變化情況下濃度計算解析式的推導，是比較困難的，只能限定侵蝕或埋藏速率的變化次數有限。埋藏速率僅一次變化的情況相對簡單又較具有現實意義，因此后文只推導埋藏速率一次變化情況下的濃度計算等式。

埋藏速率可視為負侵蝕速率，在直觀上是比較容易理解的，對于距地表x0cm深度處的樣品，侵蝕Δxcm將使其深度變為x0-Δxcm，埋藏Δxcm將使其深度變為x0+Δxcm，從而使樣品距地表深度的變化量互為相反數。

埋藏速率可視為負侵蝕速率的原理，將給各種埋藏理論模式下的濃度計算等式的推導帶來極大方便，從而無需每種埋藏理論模式最終核素濃度的計算都需要求解復雜的微分方程。

2 勻速埋藏模式下埋藏年齡的計算原理

由于采樣時樣品的當前埋藏深度是可直接測量的，所以相對于暴露過程的侵蝕量無法測量，埋藏過程相當于多了一個已知條件，因此埋藏過程相對于暴露過程，具有獨特的計算等式。下列推導的計算等式中如出現含埋藏深度“x”的參數項，則是埋藏模式特有計算等式，如不出現“x”，則不僅適用于埋藏，也適用于將侵蝕速率視為負埋藏速率的對應暴露過程。

快速深埋藏模式使用較多，此處不再討論，下面推導一些新的埋藏理論模式的濃度計算方法。

2.1 無原始濃度勻速埋藏

由于快μ介子和負μ介子對應的生成速率較散裂反應對應的生成速率隨著距地表深度增加而下降的速度要慢，因此用(8)計算勻速埋藏情況的埋藏年齡將帶來較大誤差，應對散裂反應、快μ介子和負μ介子對應的生成速率部分都應用(8)式：

注意P散、P快、P負對應的都是樣品所在深度位置而不是地表的生成速率。上式可簡化為：

(11)

可變化為包含當前樣品深度“x”而不包含埋藏年齡“T”參數的形式：

其中下標i=1對應樣品的散裂反應部分，i=2對應樣品的快μ介子部分，i=3對應樣品的負μ介子部分，可通過Pi=Pi地表e-ρx/Λi將樣品核素生成速率置換為地表核素生成速率。不同經緯度和高程的Λi值可以查閱相關文獻得到，Pi地表可通過特定經緯度和高程位置的生成速率計算方法計算獲得(Granger, 2014; Dunai, 2000, 2010; Stone, 2000; Liftonetal., 2008, 2014)。

(11)式中表面上看有兩個變量，即埋藏時間T和勻速埋藏速率β，但由于T=x/β，而埋藏深度x在采樣時可測量得到，因此變化形式后只有一個獨立變量即一個未知數β，測試單種核素濃度就可解出來勻速埋藏速率β，并進而得到埋藏年齡T。

無原始濃度勻速埋藏對應無原始濃度勻速侵蝕暴露。后者需求ε和T兩個未知數，而前者由于埋藏深度已知的條件，只需求β一個未知數，因此僅測量一種放射性核素，也可以求解出埋藏年齡。

埋藏前地表侵蝕速率極快從而樣品中的核素濃度很小，或者勻速埋藏過程持續時間長達4～5個核素半衰期因而埋藏前濃度衰變基本清零，接近無埋藏前原始濃度的假設。

2.2 有原始濃度勻速埋藏

埋藏前沒有原始濃度的可能性不是很大，考慮埋藏前存在原始濃度N0，(11)式改寫為：

(12)

變化形式得：

(12)式實質上和(1)式或(10)式是相同的，因此可視為統一的等式。可仍然將其中N0替換為假設埋藏前為穩態侵蝕(增加未知數埋藏前侵蝕速率ε)或0侵蝕速率(增加未知數埋藏前暴露時間T0)來聯立方程。在此模式下需要測試樣品的兩種放射性宇宙成因核素(如10Be和26Al)濃度，來解出埋藏速率β，以及埋藏前侵蝕速率ε(或暴露時間T0)兩個未知數，并由T=x/β求解出埋藏時間T。

快速深埋藏模式可作為此模式β非常大的特例，此時(12)式中新生濃度：

N=N0e-λT

由于計算N0需要增加一個條件(未知數)，該模式需要測量兩種放射性核素濃度。

如果是有原始濃度的快速淺埋藏，那么埋藏就位后既有繼承濃度，又有新生濃度，但侵蝕或埋藏速率變為0。

有原始濃度的勻速埋藏對應有原始濃度勻速侵蝕暴露。快速深埋藏模式對應有原始濃度而開始侵蝕速率極大，隨后侵蝕速率為0且無新生濃度的一種相對“奇特”的暴露；有原始濃度快速淺埋藏對應有原始濃度而初始侵蝕速率極大，后期侵蝕速率為0的暴露，后兩種暴露模式應屬于極不常見情況。

2.3 無原始濃度埋藏速率一次變化

在此模式下可將埋藏總時間T劃分為T1、T2兩個階段，T1階段勻速埋藏速率為β1，T2階段勻速埋藏速率為β2，此模式下將樣品的核素生成速率Pi換算為地表對應的Pi地表來計算比較方便。樣品當前濃度為第一階段生成濃度在第二階段時間衰變后的余量與第二階段的生成濃度之和。

(13)

改變形式：

此時有埋藏速率β1、β2，埋藏階段T1、埋藏階段T2四個變量需求解，即使加上x=β1T1+β2T2的限定條件，仍然需要三個獨立方程解三個未知數。因此如只能測試兩種核素，此種模式就只適用于早期勻速埋藏后期埋藏速率停止的特殊模式(即β2=0，本質上是自然埋藏速率等于自然侵蝕速率)，這時就僅需要同一樣品測試兩種核素(如10Be和26Al)來解出β1、T2兩個未知數(T1可由T1=x/β1解出)，此時等式改變為：

(14)

變化形式：

無原始濃度埋藏速率一次變化對應無原始濃度侵蝕速率一次變化暴露歷史。

2.4 有原始濃度埋藏速率變化的通用模式

該種模式應屬于自然界埋藏的最普遍模式，對應有原始濃度侵蝕速率連續變化的暴露過程，也是難以用微積分方法得到初等函數關系式的，但可以利用埋藏速率可視為負侵蝕速率的原理，以數列累加方法模擬埋藏過程中核素濃度的變化規律(黃費新等, 2020)，相應的參數ε采用負值。由于這種模式過于復雜，僅適宜以計算機編程模擬反演埋藏過程。

此種模式下可推導一些相對簡單的情況，即埋藏速率僅發生一次變化，則參考(13)式，得：

(15)

該式中有T1、T2、β1、β2，以及替換N0將產生的初始暴露時間T0或初始暴露期間侵蝕速率ε五個變量，雖然x=β1T1+β2T2可以提供一個限定條件減少方程組中的一個方程，理論上同一樣品仍需測試四種放射性核素得出四個方程方能求解，目前的實驗精度可能滿足不了求解需要。因此設第二階段埋藏速率β2為0，可以繼續簡化求解條件，從而得到：

(16)

此時有T0、T2和β1三個未知數(T1可由T1=x/β1解出)，只需測量三種放射性核素。如測試樣品為石英，可考慮測量10Be、26Al和21Ne(視為半衰期為0的放射性核素)。

如果再假設第一階段埋藏足夠快(T1=0)，就變回了有原始濃度快速淺埋藏模式：

(17)

此時就只需測量兩種核素了。

以上埋藏模式對應有原始濃度而侵蝕速率各種不同變化情況下的暴露過程。

3 誤差估計

為了簡單化，誤差估計中僅考慮濃度測試誤差引起的年齡誤差，其它參數方面的誤差將使總誤差變大。以下各式如最后埋藏年齡誤差計算中出現負值，則取絕對值得到正值。

3.1 無原始濃度勻速埋藏模式埋藏年齡的誤差

該埋藏模式的埋藏年齡由T=x/β計算得到，設x測量很準確，那么T誤差主要由N測量誤差導致的β計算誤差傳遞，即dT=-x/β2dβ。誤差dβ可根據(11)式的變化形式對β求導得到：

dβ=N(β)′dN

3.2 有原始濃度勻速埋藏模式埋藏年齡的誤差

該埋藏模式需要測定兩種核素，設為A、B，埋藏年齡仍然由T=x/β計算得到，同樣設x測量很準確，那么dT=-x/β2dβ。dβ根據(12)式的變化形式，以早期暴露期間為穩態侵蝕方案為例，組成方程組來推導：

這是有兩個變量ε、β(NA、NB是ε、β的函數)的方程組求導問題，可根據全微分公式，先求偏微分，再由偏微分得出全微分。

方程組對β求偏微分：

埋藏速率變化情況下，埋藏年齡的誤差估計的數學表達式將更為復雜，但可參考以上方法估算。或者，也可以采用數理統計的方法估計誤差，即對多個樣品測試計算結果進行相對其算術平均值的誤差統計。

4 討論

原地生成宇宙成因核素測年方法中，暴露年齡計算一般都只計算最小暴露年齡，忽略侵蝕速率的影響，在求得的最小暴露年齡在十萬年級別以上時，會產生較大的誤差，從而可能和別的測年方法得到的年齡數值無法進行對比。在埋藏年齡計算時，由于侵蝕發揮作用很小(尤其是只考慮扣除侵蝕量以后的凈埋速率時)，因而侵蝕的影響不會如計算暴露年齡時那樣對埋藏年齡計算結果產生很大的誤差。不過，由于難以確定埋藏前的準確濃度，以及難以明確埋藏的實際進行過程，即使選擇了比較合理的埋藏模式，受方法學本身的限制，計算結果可能仍然誤差較大。例如，侵蝕速率和埋藏速率實際上可能一直在持續變化，但目前的研究水平，侵蝕速率和埋藏速率的實際歷史變化情況很少得到系統性研究(有時采用模型模擬)，在暴露和埋藏年齡計算時都只能當成恒定值處理，最多假設有限次變化，而侵蝕速率和埋藏速率分別對最終暴露年齡或埋藏年齡計算結果往往具有實質性影響，有可能導致與實際情況相距甚遠。利用多樣品的埋藏年齡剖面法，有助于從數據間的對比關系判斷結果的合理性。

另外，在推導侵蝕和埋藏過程中核算濃度的計算等式中，默認了樣品對應的地表高程不變。實際上，侵蝕和埋藏都會改變樣品對應的地表高程，但一般情況下，被研究樣品所處的地表環境，樣品所在位置的侵蝕厚度往往僅在米級或十米級，而埋藏厚度則常在二十米級，再考慮地殼均衡作用的彌補，侵蝕或埋藏對樣品對應的地表高程改變量會更小，這樣小的地表高程變化對地表核素生成速率的影響非常微弱，因此對暴露或埋藏年齡計算結果影響很小，但如果侵蝕量極大或是巨厚沉積，侵蝕或埋藏過程中對應的地表核素生成速率變化將很明顯，由此將對暴露或埋藏年齡計算結果帶來較大誤差。

埋藏年齡的計算結果，顯然與采用何種埋藏模式緊密相關。究竟應該采用何種模式來計算應綜合已知資料進行分析，或者在野外調查時，對照樣品所處自然條件和周邊環境，對涉及的埋藏過程進行初步判斷，也可從多個樣品間埋藏年齡的內協性，確定采用的埋藏模式是否合理。

5 結論

(1)埋藏速率和侵蝕速率可作為互為相反數的同一參數，即埋藏速率視為負侵蝕速率，從而原地生成宇宙成因核素暴露和埋藏年齡的計算方式可以統一起來。

(2)埋藏過程除洞穴沉積等快速深埋藏模式外，應考慮埋藏是一個持續過程，因此需要考慮埋藏速率對核素最終濃度的影響，從而推導不同埋藏模式的核素濃度計算方式。

(3)埋藏速率變化下核素濃度的計算等式，也適用于對應的侵蝕速率變化下的暴露過程核素濃度的計算。