基于滲流―固體相互作用的三維水文地質(zhì)參數(shù)估計(jì)研究

［摘要］為準(zhǔn)確反映地下水系統(tǒng)的真實(shí)情況，并提高水文地質(zhì)參數(shù)的估計(jì)精度，研究提出基于三維流固耦合數(shù)學(xué)模型，對(duì)三維水文地質(zhì)參數(shù)進(jìn)行反演估計(jì)。對(duì)研究構(gòu)建的三維水文地質(zhì)參數(shù)反演算法進(jìn)行有效性研究，發(fā)現(xiàn)其導(dǎo)水系數(shù)反演誤差為5.54%，貯水系數(shù)計(jì)算誤差為2.74%，水位降深反演誤差為0.2 m，運(yùn)算時(shí)間為34.2 min，其運(yùn)算速度與精度性能較其他對(duì)比算法更高。綜上結(jié)果可知，研究構(gòu)建的三維水文地質(zhì)參數(shù)反演算法可應(yīng)用于地下水流動(dòng)特性和水文地質(zhì)參數(shù)的研究中，為未來水資源管理和土地利用規(guī)劃提供理論參考。

［關(guān)鍵詞］滲流—固體；相互作用；水文地質(zhì)；參數(shù)估計(jì)；水資源管理

地下水資源是人類生活和經(jīng)濟(jì)發(fā)展的重要組成部分。為了更好地管理和保護(hù)地下水資源，需對(duì)地下水的流動(dòng)特性和水文地質(zhì)參數(shù)進(jìn)行準(zhǔn)確估計(jì)。但地下水流動(dòng)是一個(gè)復(fù)雜得多的物理過程，其涉及水力、滲透性、彈性和化學(xué)等多個(gè)方面的相互作用[1，2]。傳統(tǒng)的水文地質(zhì)參數(shù)估計(jì)方法往往忽略了地下水系統(tǒng)中滲流和固體之間的相互作用，導(dǎo)致估計(jì)結(jié)果的不準(zhǔn)確性[3，4]。三維數(shù)學(xué)模型在描述滲流―固體相互作用中具有良好的適用性，其可基于物理方程和實(shí)驗(yàn)數(shù)據(jù)，模擬地下水的流動(dòng)和固體介質(zhì)的響應(yīng)[5]。同時(shí)通過對(duì)水文地質(zhì)參數(shù)的估計(jì)，可以更好地理解地下水流動(dòng)的機(jī)理和特性[6]。故研究提出基于固體―流體相互作用的理論，通過構(gòu)建固流耦合數(shù)學(xué)模型，并利用紐曼函數(shù)對(duì)該數(shù)學(xué)模型中的貯水率及滲透系數(shù)進(jìn)行反演，實(shí)現(xiàn)水文地質(zhì)參數(shù)的精確估計(jì)，提高數(shù)學(xué)模型的擬合能力。

1 試驗(yàn)工程概況

研究針對(duì)貴州某工程進(jìn)行現(xiàn)場(chǎng)抽水試驗(yàn)，該試驗(yàn)區(qū)域的面積約7000 m2，周長(zhǎng)為377 m。試驗(yàn)地區(qū)的底層及含水層劃分概況如圖1所示。

如圖1所示，試驗(yàn)區(qū)的地質(zhì)主要為松散巖類，地下水類型主要為松散巖類孔隙水。根據(jù)地下水的成因及年代，可將其分為淺部土層的潛水，及粉細(xì)砂中細(xì)砂層中的承壓水。實(shí)驗(yàn)場(chǎng)景內(nèi)共布置三口試驗(yàn)觀察井，首先，研究在淺部土層的潛水區(qū)域布置一口試驗(yàn)觀察井。這口試驗(yàn)井將用于觀測(cè)和采集淺部土層的潛水水位和水質(zhì)數(shù)據(jù)。隨后，研究在粉細(xì)砂中細(xì)沙層中的承壓水區(qū)域布置兩口試驗(yàn)井。這兩口試驗(yàn)井將用于觀測(cè)和采集承壓水的水位和水質(zhì)數(shù)據(jù)。承壓水是深部地下水，在巖石或土層中受到一定壓力的水。通過分析承壓水的特性和壓力變化，研究可分析地下水的流動(dòng)性和水文特征，進(jìn)而分析基于滲流―固體相互作用的三維水文地質(zhì)參數(shù)。

2 基于滲流―固體相互作用的三維水文地質(zhì)參數(shù)估計(jì)

2.1 滲流―固體相互作用的三維數(shù)學(xué)模型構(gòu)建

滲流―固體相互作用是指在多孔介質(zhì)中，流體（水或氣體）與固體顆粒之間的相互作用過程[7]。在這個(gè)過程中，流體通過多孔介質(zhì)的孔隙和通道，與固體顆粒表面發(fā)生接觸和相互作用[8]。在傳統(tǒng)的陸地工程數(shù)字仿真模型中，為保證計(jì)算的復(fù)雜度及計(jì)算時(shí)間，一般將固體假定為絕對(duì)剛體，對(duì)固流相互作用不作考慮[9]。但在地下水開采等水資源管理工程中，固體―流體的相互作用則不可被忽視。在流體―固體的相互作用中，流體和固體彼此影響，并通過相互作用來實(shí)現(xiàn)能量和質(zhì)量的傳遞，該過程又稱流固耦合[10]。此外，研究還對(duì)多孔介質(zhì)含水層中的流固耦合進(jìn)行考慮，構(gòu)建三維流固耦合數(shù)學(xué)模型，以更加準(zhǔn)確地反映地下水流動(dòng)實(shí)際運(yùn)動(dòng)情況。三維流固耦合數(shù)學(xué)模型的基本內(nèi)容如圖2所示。

如圖2，三維流固耦合數(shù)學(xué)模型包括水的壓縮性方程、多孔介質(zhì)的壓縮性方程、固體力學(xué)方程、多孔介質(zhì)的壓縮方程及多孔介質(zhì)的含水層參數(shù)。其中地下水控制方程是描述地下水彈性變形的基本方程，研究利用虎克定律和Taylor級(jí)數(shù)展開來表示其水流彈性變形，其計(jì)算公式如（1）所示。

式（1）中，ρ 表示流體密度；β 表示水體積的可變形系數(shù)；p 表示壓力。多孔介質(zhì)的含水層可能存在非彈性變形及彈性變形的情況，因此，研究利用虎克定律對(duì)其進(jìn)行計(jì)算，計(jì)算公式如（2）所示。

式（2）中，α 表示多孔介質(zhì)體積的壓縮系數(shù)；δ 表示多孔介質(zhì)骨架上的有效應(yīng)力；Vb表示多孔介質(zhì)的體積；Vs表示固體部分的體積；Vv表示空隙體積；n 表示多孔介質(zhì)的孔隙度。水和多孔介質(zhì)的壓縮性計(jì)算公式如（3）所示。

σ = λσs + （1 - λ） p （3）

式（3）中，σ 表示一維平衡狀態(tài)的垂直總壓應(yīng)力；σs表示粒間應(yīng)力的垂直分量；λ 表示單位面積中的粒間接觸面積水平投影。含水層參數(shù)包括貯水率及滲透系數(shù)，其是用于描述多孔介質(zhì)中水的分布和流動(dòng)特性的重要參數(shù)。其中貯水率是指多孔介質(zhì)中孔隙的體積與總體積之比。表示多孔介質(zhì)中水可以存儲(chǔ)的容量。貯水率越大，說明多孔介質(zhì)中可以存儲(chǔ)的水越多。貯水率的計(jì)算公式如（4）所示。

式（3）中，σ 表示一維平衡狀態(tài)的垂直總壓應(yīng)力；σs表示粒間應(yīng)力的垂直分量；λ 表示單位面積中的粒間接觸面積水平投影。含水層參數(shù)包括貯水率及滲透系數(shù)，其是用于描述多孔介質(zhì)中水的分布和流動(dòng)特性的重要參數(shù)。其中貯水率是指多孔介質(zhì)中孔隙的體積與總體積之比。表示多孔介質(zhì)中水可以存儲(chǔ)的容量。貯水率越大，說明多孔介質(zhì)中可以存儲(chǔ)的水越多。貯水率的計(jì)算公式如（4）所示。

式（4）中，μs表示貯水率；g 表示重力系數(shù)；H 表示地下水單位面積的剖面高度；dVb表示含水層體積降低排放的水量；dV 表示遇水體積增加排放的水量。滲透系數(shù)作為描述多孔介質(zhì)中水流動(dòng)性質(zhì)的參數(shù)，其可表示為單位時(shí)間內(nèi)單位面積的水通過多孔介質(zhì)的能力。滲透系數(shù)越大，說明多孔介質(zhì)中水的滲透性越好，流動(dòng)性越強(qiáng)。而導(dǎo)水系數(shù)是多孔介質(zhì)對(duì)流體的滲透能力的衡量指標(biāo)，其表示單位時(shí)間內(nèi)單位面積的介質(zhì)能夠通過的水量。導(dǎo)水系數(shù)越大，說明介質(zhì)滲透性能越好，水分能夠更快地通過。導(dǎo)水系數(shù)的大小受到介質(zhì)孔隙度、孔隙連通性、孔隙大小和土壤顆粒組成等因素的影響。研究利用滲透試驗(yàn)對(duì)導(dǎo)水系數(shù)的值進(jìn)行確定。滲透系數(shù)和導(dǎo)水系數(shù)的計(jì)算公式如式（5）所示。

T = KM （5）

式（5）中，T 表示導(dǎo)水系數(shù)；M 表示滲透系數(shù)。隨后，研究基于以上參數(shù)進(jìn)行固流耦合的數(shù)學(xué)模型的構(gòu)建。

2.2 基于流固耦合模型的三維水文地質(zhì)參數(shù)反演

完成固流耦合數(shù)學(xué)模型的構(gòu)建后，研究利用紐曼函數(shù)對(duì)數(shù)學(xué)模型中的貯水率及滲透系數(shù)進(jìn)行反演，獲取三維水文地質(zhì)參數(shù)，以更加精確地描述多孔介質(zhì)中水的流動(dòng)變化。紐曼函數(shù)標(biāo)準(zhǔn)曲線分為A組曲線及B組曲線，分別代表標(biāo)準(zhǔn)曲線的上凹部分及下凸部分[11，12]。其反演的基本思想為通過對(duì)A組和B組的曲線與觀察井的實(shí)際時(shí)間―降深曲線（s-t曲線）進(jìn)行最優(yōu)配線，以獲取最優(yōu)水文地質(zhì)變量[13]。但由于固流耦合模型的復(fù)雜性，傳統(tǒng)的優(yōu)化算法往往面臨著維度高、非線性強(qiáng)等問題，導(dǎo)致求解難度較大。而隨機(jī)權(quán)重粒子群優(yōu)化算法（Random Weighted Particle Swarm Optimization，RandWPSO）采用了粒子群優(yōu)化的思想，通過模擬鳥群的行為來搜索最優(yōu)解，能夠有效地應(yīng)對(duì)流固耦合模型的高維度和非線性問題，提高優(yōu)化的效率[14，15]。故研究選擇RandWPSO算法作為流固耦合模型的優(yōu)化配線方法，其配線流程如圖3所示。

如圖3所示，研究首先對(duì)RandWPSO算法的粒子維度、粒子種群大小、隨機(jī)權(quán)重平均值的范圍及最大迭代次數(shù)等參數(shù)進(jìn)行設(shè)置。隨后，研究對(duì)粒子群的速度及位置進(jìn)行初始化操作，并利用適應(yīng)度函數(shù)對(duì)粒子群進(jìn)行評(píng)價(jià)并對(duì)比，獲取最小適應(yīng)度的粒子作為群體最優(yōu)粒子，以此更新粒子速度及位置，直至到達(dá)最大迭代次數(shù)。隨后，研究對(duì)B組進(jìn)行配線，獲取最優(yōu)曲線的配線參變量。紐曼函數(shù)的配線方式包括兩種方式，研究分為方法1及方法2。以B組曲線為例，方法1的具體過程表示為，首先對(duì)s-t曲線進(jìn)行繪制，并將該曲線的后半段對(duì)B組曲線進(jìn)行擬合，并確定其最優(yōu)配線參變量。隨后選取任意擬合點(diǎn)，計(jì)算B組導(dǎo)水系數(shù)和導(dǎo)水度。配線參變量通過計(jì)算曲線的導(dǎo)數(shù)和二階導(dǎo)數(shù)，找到曲線上的極值點(diǎn)和拐點(diǎn)來實(shí)現(xiàn)。根據(jù)最優(yōu)配線參變量，即可獲得與B組曲線相反的曲線，即A組曲線。將A組曲線與s-t曲線的前半段進(jìn)行最優(yōu)擬合，計(jì)算A組導(dǎo)水系數(shù)和貯水系數(shù)。方法2與方法1除先后順序外，其計(jì)算方法均一致。方法2首先對(duì)A組曲線進(jìn)行配線，獲取其配線參變量，隨后以該參變量進(jìn)行B組的配線。A組和B組的導(dǎo)水系數(shù)應(yīng)大體相近。研究將對(duì)方法1和方法2獲取的AB兩組導(dǎo)水系數(shù)差值進(jìn)行計(jì)算，選取性能更優(yōu)的配線方法作為參數(shù)的計(jì)算方法。最后，研究通過A組和B組配線的曲線及結(jié)果圖，輸出水文地質(zhì)參數(shù)值。

3 基于三維流固耦合數(shù)學(xué)模型的水文地質(zhì)參數(shù)反演算法實(shí)證分析

3.1 三維流固耦合數(shù)學(xué)模型有效性驗(yàn)證

為驗(yàn)證研究構(gòu)建的三維固耦合模型的有效性，研究利用抽水試驗(yàn)中觀察井的降深―時(shí)間觀測(cè)值作為數(shù)據(jù)集，對(duì)三維流固耦合模型進(jìn)行有效性驗(yàn)證。對(duì)比模型為一維流固耦合數(shù)學(xué)模型和二維流固耦合數(shù)學(xué)模型。性能評(píng)價(jià)指標(biāo)為與實(shí)際水位降深值的偏差及水位降深值的計(jì)算精度。實(shí)驗(yàn)環(huán)境為多物理場(chǎng)仿真軟件COMSOL。該三維流固耦合數(shù)學(xué)模型的水位降深模擬數(shù)據(jù)及實(shí)測(cè)水位降升對(duì)比結(jié)果如圖4所示。

圖4（a）為各三維固流耦合模型的擬合準(zhǔn)確度對(duì)比結(jié)果。如圖4（a）所示，研究提出的三維固流耦合數(shù)學(xué)模型計(jì)算出的水位降深值與實(shí)際值更為接近，其與實(shí)際值的最大偏差為0.06 m，較一維流固耦合模型的最大偏差低0.09 m，較二維固流模型的最大偏差低0.07 m，故研究提出的三維固流模型擬合準(zhǔn)確度更優(yōu)。圖4（b）為各三維固流耦合模型的擬合精確度對(duì)比結(jié)果。如圖4（b）所示，研究提出的基于固流耦合模型計(jì)算出的水位降深值計(jì)算精確度為90.1%，較其他模型更高。綜上結(jié)果可知，研究提出的固流耦合模型擬合準(zhǔn)確度及精度更高。

3.2 水文地質(zhì)參數(shù)反演算法有效性驗(yàn)證

為選擇參數(shù)反演算法中，性能更優(yōu)的配線方法。研究分別對(duì)配線方法1及配線方法2進(jìn)行計(jì)算結(jié)果分析。實(shí)驗(yàn)數(shù)據(jù)集為抽水試驗(yàn)的觀測(cè)值，實(shí)驗(yàn)環(huán)境為MATLABE。兩種方法的計(jì)算結(jié)果如表1所示。

如表1所示，方法1中，A組的導(dǎo)水系數(shù)為176.73，B組的導(dǎo)水系數(shù)為130.82，A組與B組的導(dǎo)水系數(shù)差值為45.91。在方法2中，A組的導(dǎo)水系數(shù)為176.28，B組的導(dǎo)水系數(shù)為137.55，A組與B組的導(dǎo)水系數(shù)差值為39.73。綜上結(jié)果可知，使用配線方法2的導(dǎo)水系數(shù)差值較配線方法1的導(dǎo)水系數(shù)低，故配線方法2更符合實(shí)際情況，可應(yīng)用于水文地質(zhì)參數(shù)的反演中。最后，為驗(yàn)證研究提出的基于固流耦合的參數(shù)反演算法性能，研究基于COMSOL及MATLABE對(duì)抽水試驗(yàn)觀察井的觀測(cè)值進(jìn)行計(jì)算。對(duì)比方法為基于一維固流耦合的參數(shù)反演算法、基于二維流固耦合的參數(shù)反演算法及基于三維非流固耦合的參數(shù)反演法。對(duì)比指標(biāo)為導(dǎo)水系數(shù)、貯水系數(shù)、觀察井的水位降深及運(yùn)行時(shí)間。各反演算法的計(jì)算結(jié)果對(duì)比如圖5所示。

圖5（a）為各反演算法的擬合準(zhǔn)確度對(duì)比結(jié)果。如圖5（a）所示，研究提出的基于三維耦合模型的參數(shù)反演算法導(dǎo)水系數(shù)反演誤差最小，為5.54%，較其他反演算法的誤差更低。此外，該算法的貯水系數(shù)計(jì)算誤差為2.74%，其性能較其他模型更優(yōu)。圖5（b）為各反演算法的水位降深與運(yùn)行效率對(duì)比結(jié)果。如圖5（b）所示，研究提出的基于三維耦合模型的參數(shù)反演算法水位降深反演誤差最小，為0.2 m，較其他模型的誤差更小，此外，該算法的運(yùn)算時(shí)間為34.2 min，較一維和二維的反演算法較高，但較三維非固流耦合低52 min。綜上結(jié)果可知，研究構(gòu)建的基于三維固流耦合的參數(shù)反演算法可平衡水文地質(zhì)參數(shù)的求解性能及求解速度，可對(duì)實(shí)際水位進(jìn)行預(yù)測(cè)和求解。

4 結(jié)論

為更加精確地模擬地下水流動(dòng)和貯水特性，研究利用滲流―固體相互作用構(gòu)建三維流固耦合數(shù)學(xué)模型，并將該模型與RandWPSO算法進(jìn)行融合，構(gòu)建基于滲流―固體相互作用的三維水文地質(zhì)參數(shù)反演算法，獲取水文地質(zhì)參數(shù)。對(duì)研究構(gòu)建的三維流固耦合數(shù)學(xué)模型進(jìn)行有效性驗(yàn)證，發(fā)現(xiàn)水位降深誤差為0.06 m，精確度為90.1%，模擬準(zhǔn)確度及精確度性能較其他模型更優(yōu)。此外，研究還對(duì)基于三維耦合模型的參數(shù)反演算法進(jìn)行有效性驗(yàn)證，發(fā)現(xiàn)其導(dǎo)水系數(shù)反演誤差為5.54%，貯水系數(shù)計(jì)算誤差為2.74%，水位降深反演誤差為0.2 m，顯著低于其他對(duì)比算法。此外，該算法的運(yùn)算時(shí)間為34.2 min，較三維非固流耦合低52 min。綜合以上結(jié)果可知，研究提出的基于三維耦合模型的參數(shù)反演算法可對(duì)三維水文地質(zhì)參數(shù)進(jìn)行準(zhǔn)確估計(jì)，并為地下水流動(dòng)特性和水文地質(zhì)參數(shù)的研究提供理論參考，進(jìn)而對(duì)地下水資源管理和環(huán)境保護(hù)等領(lǐng)域提供有益的支持。

[參考文獻(xiàn)]

［1］高成路，李術(shù)才，周宗青，等. 巖體水力壓裂應(yīng)力―滲流耦合近場(chǎng)動(dòng)力學(xué)模擬［J］. 同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2022，50（4）：470-480.

［2］張培森，趙成業(yè)，侯季群，等. 溫度―應(yīng)力―滲流耦合條件下紅砂巖滲流特性試驗(yàn)研究［J］. 巖石力學(xué)與工程學(xué)報(bào)，2020，39（10）：1957-1974.

［3］樊輝，呂義清，鄭劭鐸. 基于滲流―應(yīng)力耦合的采動(dòng)堆積層滑坡變形分析［J］. 煤礦安全，2020，51（11）：281-289.

［4］劉慶生，肖浩，談成亮，等. 離子型稀土礦浸出過程的滲流―反應(yīng)―應(yīng)力耦合模擬研究［J］. 中國(guó)稀土學(xué)報(bào)，2022，40（5）：880-892.

［5］慎國(guó)強(qiáng)，劉浩杰，王振濤，等. 致密儲(chǔ)層復(fù)雜縫網(wǎng)流固耦合滲流規(guī)律及產(chǎn)能模擬［J］. 大慶石油地質(zhì)與開發(fā)，2023，42（4）：157-165.

［6］周亞醒，陳浩，張飛，等. 群孔抽水試驗(yàn)方法求解含水層水文地質(zhì)參數(shù)［J］. 水文，2022，42（6）：78-82.

［7］劉淑惠，肖長(zhǎng)來，梁秀娟. 基于線性、非線性規(guī)劃的雙遼市水文地質(zhì)參數(shù)計(jì)算方法［J］. 水電能源科學(xué)，2020，38（4）：72-75.

［8］Brunetti C，Bianchi M，Pirot G，Linde N. Hydrogeological Model Selection Among Complex Spatial Priors［J］. Water Resources Research，2019，55（8）：6729-6753.

［9］樊冬冬，劉天任，譚勇，等. 富水砂層水文地質(zhì)參數(shù)反演及基坑降水圍護(hù)優(yōu)化［J］. 建筑技術(shù)，2023，54（13）：1587-1591.

［10］王旭東，呂樂，時(shí)俊，等. 抽水試驗(yàn)反演水文地質(zhì)參數(shù)的多線全程加權(quán)法［J］. 水利學(xué)報(bào)，2020，51（3）：276-285.

［11］張淵博，鄒德旋，張春韻，等. 自適應(yīng)慣性權(quán)重的粒子群優(yōu)化算法［J］. 計(jì)算機(jī)仿真，2023，40（4）：350-357.

［12］Marshall S K，Cook P G，Simmons C T，Konikow L F，Dogramaci S. An Approach to Include Hydrogeologic Barriers With Unknown Geometric Properties in Groundwater Model Inversions［J］. Water Resources Research，2022，58（7）：1-22.

［13］石鴻蕾，郝奇琛，邵景力，等. 基于多源數(shù)據(jù)的弱透水層水文地質(zhì)參數(shù)反演研究——以呼和浩特盆地某淤泥層為例［J］. 水文地質(zhì)工程地質(zhì)，2021，48（2）：1-7.

［14］李磊，陳干，唐沛，等. 基于解析法和數(shù)值法反演哈爾濱漫灘區(qū)水文地質(zhì)參數(shù)［J］. 城市軌道交通研究，2021，24（10）：54-59.

［15］邱淑偉，吳亞敏，柯昱琪，等. 基于遍歷搜索算法的水文地質(zhì)參數(shù)優(yōu)化求解［J］. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2020，50（6）：1854-1861.