考慮毛細(xì)飽水帶的非承壓含水層中仿Theis井流模型及解析解

摘要：

受制于經(jīng)典飽和滲流理論中自由面定位在潛水面處的局限性，非承壓含水層中仿Theis井流模型存在理論基礎(chǔ)不堅實、參數(shù)物理意義不明確等問題。為克服上述局限性，首先將自由面由潛水面移動到進氣值面，以定位在進氣值面處的Boussinesq方程為控制方程，修正仿Theis井流模型，并推導(dǎo)模型的解析解；之后，構(gòu)建參數(shù)反演模型，以校準(zhǔn)模型參數(shù)，并通過分別將修正模型計算的降深和經(jīng)典模型計算的理論降深與實測降深進行對比，檢驗修正模型的合理性；然后，利用參數(shù)敏感性分析探討修正模型下進氣值面高程的變化特征；最后，討論進氣值面處鉛直方向滲流分速度的計算方法和變化特征。結(jié)果表明：修正模型計算的降深曲線與實測降深曲線基本吻合；進氣值面高程隨與井中心的距離、給水度增大而增大，隨抽水流量增大而減小。在抽水早期，進氣值面高程隨滲透系數(shù)的增大而減小，在抽水后期則相反。基于進氣值面處水量均衡關(guān)系（基于線性或非線性方程，采用完整解或近似解）及進氣值面處滲流連續(xù)性方程（基于線性方程，采用完整解）等5種情況，推導(dǎo)的鉛直方向滲流分速度的解析表達計算結(jié)果基本吻合，曲線趨勢一致，表現(xiàn)為隨抽水時間的延長呈非線性增長。

關(guān)鍵詞：

進氣值面；井流模型；完整井；Boussinesq方程；解析解；參數(shù)反演模型

doi：10.13278/j.cnki.jjuese.20230151

中圖分類號：P641.2

文獻標(biāo)志碼：A

Modified Theis Well Model and Its Analytical Solution in Unconfined Aquifers Considering Capillary Saturated Zone

Jin Haoying^1，2，Cheng Dawei^"1，2，Zhan Hongbing³，Yang Shengke^1，2，Zhang Xiaofei⁴，Zhang Lin^1，2

1. School of Water and Environment， Chang’an University， Xi’an 710064，China

2. Key Laboratory of Subsurface Hydrology and Ecological Effects in Arid Region （Chang’an University），Ministry of Education，"Xi’an 710064，China

3. Department of Geology and Geophysics， Texas Aamp;M University， College Station， TX 77843-3115， USA

4. Xi’an Qujiang Meibei Lake Investment and Construction Co.， Ltd.， Xi’an 710300，China

Abstract：

Due to the limitation of setting the upper boundary at the free surface located at the water table， the Theis well model in unconfined aquifers has some problems such as unsound theoretical basis and unreasonable physical meaning of parameters. To overcome these issues， firstly， the well model is revised by moving the free surface from the water table to the air entry plane， and takes the Boussinesq equation located at the air entry plane as the governing equation， and derives the analytical solution of the model; Then， a parameter inversion model is constructed to calibrate the model parameters. By comparing the revised model drawdown and the theoretical drawdown of the classical model with the measured drawdown， it verifies the rationality of the well revised model proposed in this paper. Then， the variation characteristics of the air entry plane elevation under the modified well model are discussed by using parameter sensitivity analysis. Finally， the computational method and variation characteristics of seepage velocity in the vertical direction of the air entry plane are discussed. It is found that the drawdown curve obtained by the revised model coincides with the measured drawdown curve; The air entry plane elevation increases with increasing distance from the center of the well and the specific yield and decreases with increasing pumping flow rate. In the early stage of pumping， the air entry plane elevation decreases with increasing hydraulic conductivity， but the opposite is true in the late stage of pumping. The calculation results of the analytical expressions of the seepage velocity in the vertical direction based on the water balance relation at the air entry plane （based on the linearized or non-linearized equation， using complete solutionor approximate solution） and based on continuity equation of seepage flow at the air entry plane （based on the linearized equation， using complete solution） under five conditions are basically consisent， and the trends of the curve are consistent with each other， showing a nonlinear increase with the prolongation of pumping time.

Key words：

air entry plane；well model；fully penetrating well；Boussinesq equation；analytical solution；parameter inversion model

0"引言

非穩(wěn)定井流理論能夠揭示抽水或注水?dāng)_動下地下水水位的變化及其運動的基本規(guī)律^[1]。構(gòu)建非穩(wěn)定井流模擬地下水流動態(tài)，評估地下水水量變化，并評價其影響，對地下水資源開采、地下水環(huán)境保護具有重要意義^[2]。此外，以非穩(wěn)定井流理論為基礎(chǔ)的抽水試驗及參數(shù)反演模型是獲取含水層水文地質(zhì)參數(shù)的重要手段之一，在地下水流模型的構(gòu)建中發(fā)揮著重要作用^[3-6]。

在地下水井流問題的研究中，Theis^[7]提出的針對承壓含水層的非穩(wěn)定井流模型（Theis解）是經(jīng)典解答之一，也是井流理論的基礎(chǔ)。后繼學(xué)者在Theis井流模型的基礎(chǔ)上，針對不同的水文地質(zhì)條件，從介質(zhì)的異質(zhì)性、井儲效應(yīng)、水與介質(zhì)骨架的壓縮性等方面對Theis井流模型進行了改進^[8-13]。當(dāng)把Theis井流模型應(yīng)用于非承壓含水層時，常被稱為仿Theis井流，其區(qū)別是用給水度替換釋水系數(shù)^[14]。仿Theis井流公式被廣泛應(yīng)用到水文地質(zhì)、工程地質(zhì)中進行參數(shù)估計和預(yù)測水位降深^[15-20]。此外，學(xué)者們依據(jù)對井流過程機理的認(rèn)識，從滯后排水、介質(zhì)異質(zhì)性、補給條件等方面改進了仿Theis井流模型^[21-23]。在描述非承壓含水層井流時，仿Theis井流模型刻畫的是均質(zhì)各向同性含水層內(nèi)，完整單井定流量抽水所引起的地下水滲流動態(tài)過程^[14]。仿Theis井流模型因循經(jīng)典飽和滲流理論的相關(guān)結(jié)論，將含水層的上邊界（自由面）定位在潛水面處，該模型的控制方程為自由面定位在潛水面處的Boussinesq方程。

進氣值面是指支持毛細(xì)帶的上邊界，高于潛水面，是飽和含水量與非飽和含水量的分界面。Bear^[24]認(rèn)為自由面通常指潛水面，但從飽和-非飽和邊界的角度來考察，也可指進氣值面。然而，Cheng等^[25]論證了基于潛水面定義的潛水面方程是不自洽的。主要表現(xiàn)為：1）基于突變界面假定的潛水面方程，在假定中忽略了潛水面之上毛細(xì)飽和帶和毛細(xì)非飽和帶的影響，不符合物理真實。2）基于修正給水度的潛水面方程，潛水面處水量均衡關(guān)系不匹配。3）潛水面方程的一般式，當(dāng)考慮潛水面之上存在毛細(xì)飽和帶時，其給水度值恒為0。在滿足Dupuit假定的情況下，Boussinesq方程可視為非承壓含水層飽和滲流控制方程和潛水面方程的簡化形式，受自由面定位在潛水面處的局限性制約，具有明顯的局限性。為了克服定位在潛水面處Boussinesq方程的局限性，程大偉^[26]建立了定位在進氣值面處的Boussinesq方程。由于定位在進氣值面處的自由面使飽和滲流區(qū)域同時包含了毛細(xì)飽和帶和重力飽和帶，因而定位在進氣值面處的Boussinesq方程具備了描述毛細(xì)飽和帶滲流的功能。就仿Theis井流問題而言，這實質(zhì)上賦予仿Theis井流模型新的上邊界和控制方程。

本文首先對仿Theis井流模型進行修正，將自由面由潛水面移動到進氣值面，采用定位在進氣值面處的Boussinesq方程作為控制方程，并通過線性變換及拉普拉斯變換，推導(dǎo)修正仿Theis井流模型的解析解；其次，基于新的解析解構(gòu)建參數(shù)反演模型，以校準(zhǔn)模型參數(shù)；之后，分別將修正模型計算的降深和經(jīng)典模型計算的理論降深與實測降深進行對比，檢驗修正后模型的合理性；然后，利用參數(shù)敏感性分析探討修正模型下進氣值面高程的變化特征；最后，討論了進氣值面處鉛直方向的滲流分速度的計算方法和變化特征。旨在探索考慮毛細(xì)飽水帶的仿Theis井流問題的解答，以豐富和發(fā)展井流相關(guān)理論。

1"滲流模型

1.1"數(shù)學(xué)模型

非承壓含水層的自由面應(yīng)定位在進氣值面處，該處的壓力水頭為進氣值水頭^[25]，相應(yīng)的自由面定位在進氣值面處的完整井流模型如圖1所示。當(dāng)采用定位在進氣值面的Boussinesq方程來刻畫該滲流過程時，需要對經(jīng)典仿Theis井流模型的基本假定做如下修改與補充：飽和滲流區(qū)的上邊界在進氣值面處，初始時刻進氣值面是水平的，稱為靜進氣值面，抽水時降低的進氣值面稱為動進氣值面；進氣值面降深等于靜進氣值面與動進氣值面之差，記作s_ha（r，t），如圖1所示。

程大偉^[26]建立的定位在進氣值面處的Boussinesq方程如式（1）所示：

x（K_xhhx）+y（K_yhhy）+W=μ_haht。（1）

式中：x、y為直角坐標(biāo)系下的坐標(biāo)值，m；K_x、K_y分別為x和y方向的飽和滲透系數(shù)，m/min；h為含水層厚度（由隔水底板到進氣值面的距離），m；W為進氣值面處的入滲補給率，m/min；μ_ha為進氣值面處的給水度，定義為在橫截面為單位面積的巖土柱內(nèi)，在重力和毛細(xì)力的共同作用下，進氣值面移動的單位高度在進氣值面移動范圍內(nèi)水量的變化量^[25]；t為抽水時間，min。

自由面定位在進氣值面完整井的數(shù)學(xué)模型如下：

²ur²+1rur=μ_haTut "（rgt;0，tgt;0） ；""（2a）h（r，0）=h₀"""（t=0，rgt;0） ；"""""""（2b）h（，t）=h₀""""""""""（tgt;0，r→） ；""""""（2c）limr→0（rhr）=0 "（tgt;0，h≥zgt;h+ha，r→0） ；Q2πKh_wt"（tgt;0，z≤h+ha，r→0）。

（2d）

式中：u為非線性向線性轉(zhuǎn)變的參數(shù)，u=h²/2；T為導(dǎo)水系數(shù)，m²/min，T=Kh_cp，h_cp為含水層平均厚度，在上述假定條件下，可近似取h₀；h₀為初始非承壓含水層厚度，即初始時刻從底板到進氣值面的高度，其值是初始潛水面處的高程h_0，wt與進氣值水頭h_a之差，m；Q為完整井的抽水流量，m³/min；K為飽和滲透系數(shù)，m/min；h_wt為隔水底板到潛水面的高度，m。

式（2a）為式（1）經(jīng)線性化處理并忽略進氣值面處的入滲補給率W，在極坐標(biāo)與軸對稱滲流條件下的形式，且含水層是均質(zhì)各向同性的，定位在進氣值面處的Boussinesq方程。

由于進氣值面到潛水面的距離相對于井壁處正壓區(qū)的范圍要小得多，為了簡化求解，將抽水量平均到整個飽和區(qū)（包含正壓區(qū)和進氣值面到潛水面處的負(fù)壓區(qū)），則式（2d）可寫作

limr→0（rhr）=Q2πKh（tgt;0，z≤h+ha，r→0）。 （3）

若令降深為s_ha=h₀－h(huán)，修正降深為

s′_ha=s_ha－s_ha²2h₀，則

u=h²2=12（h₀－s_ha）²=h₀²2－h(huán)₀s′_ha。（4）

在極坐標(biāo)與軸對稱滲流條件下，分別對控制方程和井流量邊界條件進行變換，可建立以s′_ha為因變量的定位在進氣值面處非承壓含水層的完整井?dāng)?shù)學(xué)模型，其中井壁處的邊界條件對應(yīng)于式（3）：

T²s′_har²+1rs′_har=

μ_has′_hat（tgt;0，rgt;0） ；s′_ha（r，0）=0 """"（t=0，rgt;0） ；s′_ha（，t）=0 ""（tgt;0，r→） ；limr→0（rs′_har）=－Q2πT """（tgt;0，r→0） 。（5a）

（5b）

（5c）

（5d）

1.2"模型求解

為了便于求解數(shù)學(xué)模型，可對式（5）作關(guān)于t的Laplace變換，得

d²s′_hadr²+1rds′_hadr=Δas′_ha"； """"""s′_ha（，Δ）=0 ； """""""""""""""""""""""""limr→0rs′_ha（r，Δ）r=－Q2πTΔ 。 "（6a）

（6b）

（6c）

式中：s′_ha為s′_ha關(guān)于t的Laplace變換的象函數(shù)；Δ為Laplace算符；a為導(dǎo)壓系數(shù)，a=T/μ_ha，

m²/min。

通過作變量置換及拉普拉斯逆變換，可得修正降深的表達式為

s′_ha（r，t）=Q4πTW（u_μ）。（7）

式中：u_μ=r²μ_ha4Tt；W（u_μ）=∫u_μ1ye^－ydy（u_μgt;0）。

若換成實際降深公式，則實際降深的完整解為

s_ha（r，t）=h₀－h(huán)₀²－Q2πKW（u_μ）。 （8）

進氣值面高程的完整解為

h（r，t）=h₀²－Q2πKW（u_μ）。（9）

當(dāng)u_μ=0.01～0.05時，有

W（u_μ）≈ln2.25Ttr²μ_ha。 （10）

因此，將式（10）代入式（8），可得進氣值面處水位降深的近似解為

s_ha（r，t）=h₀－h(huán)₀²－Q2πKln2.25Ttr²μ_ha。（11）

相應(yīng)地，進氣值面高程的近似解為

h（r，t）=h₀²－Q2πKln2.25Ttr²μ_ha。 （12）

潛水面處水位高程h_q的近似解為

hq（r，t）=ha+h₀²－Q2πKln2.25Ttr²μ_ha。（13）

潛水面處水位降深s_hq的近似解為

s_hq（r，t）=h₀－h(huán)₀²－Q2πKln2.25Ttr²μ_ha。（14）

1.3"構(gòu)建參數(shù)反演模型

在上述解析解的基礎(chǔ)上構(gòu)建相應(yīng)的參數(shù)反演模型，可以用于進氣值面處的給水度μ_ha、飽和滲透系數(shù)K、進氣值水頭h_a、含水層的平均厚度h_cp、初始潛水面處的高程h_0，wt等參數(shù)的校準(zhǔn)。反演模型包括目標(biāo)函數(shù)和約束條件。

目標(biāo)函數(shù)為

minf=∑i∑j（s_hq，cij－s_hq，tij）²。 "（15a）

約束條件為

μ_ha，minlt;μ_halt;μ_ha，max；（15b）

K_minlt;Klt;K_max；（15c）

ha_，minlt;halt;h_a，max；（15d）

h_cp，minlt;h_cplt;h_cp，max； （15e）

h_{0，wt，min}lt;h_0，wtlt;h_{0，wt，max}。（15f）

式中：s_hq，cij和s_{hq，tij}分別為第j個時刻和第i個觀測點的水位降深的計算值和觀測值，m，其中s_hq，tij可通過式（14）來獲得；min與max表示反演參數(shù)的下限與上限；μ_ha的上限為飽和含水量與殘余含水量之差，下限為0；K的上、下限根據(jù)沉積物類型確定，m/min；h_a的上、下限根據(jù)沉積物類型確定，m；h_cp和h_0，wt的上下限可由其初始值的±40%范圍決定，m。

2"模型驗證

為了驗證修正后仿Theis井流模型解析解的合理性，將本文提出的解析解計算得到的潛水面處降深（式（14））與前人實測的潛水面處的降深進行對比。江峰等^[27]在福建某濱海風(fēng)積-海積平原進行了10組抽水試驗，10組試驗分別對應(yīng)10口不同位置的抽水井，本文選取了其中三口不同抽水井抽水試驗數(shù)據(jù)，標(biāo)記為a、b、c井（與井中心的距離r分別為2.60、1.90、3.50 m），用作本研究的實測數(shù)據(jù)。各組抽水試驗的基本情況如表1所示，其中滲透系數(shù)、潛水面處的給水度等是利用Aquifer Test 3.0軟件基于經(jīng)典仿Theis井流模型開展參數(shù)反演獲取的估計值^[27]，稱為經(jīng)典模型反演值；潛水面處實測降深數(shù)據(jù)如圖2所示。基于此實測數(shù)據(jù)，利用1.3節(jié)構(gòu)建的參數(shù)反演模型對參數(shù)進行反演。待反演的參數(shù)及其初始值與上下限見表2，反演結(jié)果如表3所示，稱為修正模型反演值。將校準(zhǔn)后的參數(shù)值（表3）代入式（14）中，可獲得基于修正模型計算的潛水面降深，在理論降深計算中使用表1中的參數(shù)，結(jié)果見圖2。由圖2可知，理論降深與實測降深的差異很大。

根據(jù)江峰等^[27]的實測降深數(shù)據(jù)，利用式（15）對修正的仿Theis模型（式（14））進行參數(shù)校準(zhǔn)，獲得的結(jié)果如表3所示。由表3可知，在a、b、c三井的初始潛水面高程反演值中，a井相對誤差最小，為0.2%；b井相對誤差最大，為24.4%。在a、b、c三井的滲透系數(shù)反演值中，b井的相對誤差最大，為273.3%；c井相對誤差最小，為84.6%。通過對比修正模型計算降深和實測降深（圖2）可知，修正模型計算降深與實測降深基本吻合，且隨著抽水時間的延長，修正模型計算降深曲線與實測降深曲線幾乎重合。原因在于水位降深變化速率會隨時間的延長而逐漸減小，降深曲線接近平緩，此時抽水過程更符合仿Theis井流模型的假定條件。通過上述對比還表明本論文所提出的解析解是合理的。

將表3獲得的參數(shù)反演結(jié)果代入式（12），計算得到a、b、c三井抽水試驗對應(yīng)的進氣值面高程，如圖3所示。由圖3可見，當(dāng)r分別固定（a、b、c井，r=2.60、1.90、3.50 m）時，隨抽水時間的延長，r處進氣值面高程呈非線性降低的趨勢。在抽水試驗前期，進氣值面高程急速下降；隨抽水時長增加，進氣值面高程下降速率減小，進氣值面高程曲線趨于平緩。

3"討論

3.1"參數(shù)敏感性

修正的仿Theis井流模型的進氣值面高程表達式為式（12），在式（10）的條件下，將式（12）對r、Q、μ_ha、K分別求偏導(dǎo)，通過參數(shù)敏感性分析來探討各水文地質(zhì)參數(shù)對完整井在進氣值面高程的影響。

對修正的仿Theis井流模型在進氣值面高程h（r，t）取r的一階導(dǎo)數(shù)及二階導(dǎo)數(shù)，有

h（r，t）r=Q2πKr（h₀²－Q2πKln2.25Ttr²μ_ha）^－12gt;0； （16）

²h（r，t）r²=－Q2πKr²（h₀²－Q2πKln2.25Ttr²μ_ha）^－12－Q²4π²K²r²（h₀²－Q2πKln2.25Ttr²μ_ha）^－32lt;0。（17）

由式（16）（17）可知，對于給定的t時刻r處的h（r，t）與r的關(guān)系曲線在剖面上是一條單調(diào)遞增凸曲線，且當(dāng)r越小，曲線在r方向上坡度的絕對值越大；r越大，坡度的絕對值越小；當(dāng)r→時，h（r，t）r→0，說明在無窮遠(yuǎn)處進氣值面高程曲線與初始水頭線相切。

對h（r，t）取Q的一階導(dǎo)數(shù)，有

h（r，t）Q=

－14πK（h₀²－Q2πKln2.25Ttr²μ_ha）^－12ln2.25Ttr²μ_halt;0。（18）

此時，在其他水文地質(zhì)參數(shù)不變的條件下，給定t時刻r處的h（r，t）隨Q的增大而減小。

對h（r，t）取μ_ha的一階導(dǎo)數(shù)，有

h（r，t）μ_ha=Q4πKμ_ha（h₀²－Q2πKln2.25Ttr²μ_ha）^－12gt;0。 （19）

此時，h（r，t）與μ_ha成正比。在其他水文地質(zhì)參數(shù)不變的條件下，給定t時刻r處的

h（r，t）隨μ_ha

的增大而增加。

對h（r，t）取K的一階導(dǎo)數(shù)，有

h（r，t）K=Q4πK²（h₀²－Q2πKln2.25Ttr²μ_ha）^－12（ln2.25Kh₀tr²μ_ha－1）。

（20）

當(dāng)Kgt;er²μ_ha2.25h₀t時，表明h（r，t）與K成正比，反之，h（r，t）與K成反比。由于er²μ_ha2.25h₀t是關(guān)于時間t的反比例函數(shù)（隨時間的增大，該函數(shù)值由+向0遞減），這意味著總存在某個時刻t_K，使得K=er²μ_ha2.25h₀t_K。當(dāng)tlt;t_K時，Klt;er²μ_ha2.25h₀t_K，結(jié)合式（20）可知，h（r，t）與K成反比；當(dāng)tgt;t_K時，Kgt;er²μ_ha2.25h₀t_K，可知h（r，t）與K成正比。綜合上述分析可知：在其他水文地質(zhì)參數(shù)不變的條件下，在抽水早期，給定t時刻r處的h（r，t）隨K的增大而減小；在抽水后期，給定的t時刻r處的h（r，t）隨K的增大而增大。

3.2關(guān)于進氣值面處鉛直方向滲流分速度的討論

學(xué)界通常認(rèn)為基于Boussinesq方程無法獲取鉛直方向的滲流分速度^[28]，這種觀點是值得商榷的。原因在于：Dupuit假定表述為當(dāng)滲流的鉛直分速度與水平分速度相比顯得很小時，可用自由面處的滲流水平分速度近似代替整個鉛直線上的水平分速度^[14]，該假定并未否認(rèn)鉛直方向滲流分速度的存在。在Boussinesq方程的推導(dǎo)過程中承認(rèn)存在鉛直方向滲流分速度，并基于滲流連續(xù)性方程和水量均衡關(guān)系分別獲取鉛直方向滲流分速度v_z，將二者聯(lián)立獲得v_z的隱式表達，得到了Boussinesq方程^[14]。因此，基于Boussinesq方程可以求得v_z。

定位在進氣值面處的Boussinesq方程推導(dǎo)的基本思路與前人^[14]相似。二者的主要區(qū)別是經(jīng)典Boussinesq方程定位在潛水面處。因此，基于定位在進氣值面處的Boussinesq方程也可以求得v_z。以下我們將進氣值面高程的完整解和近似解分別代入線性（式（5））或非線性（式（1））的定位在進氣值面處Boussinesq方程等號的左右兩側(cè)（等號左側(cè)代表基于進氣值面處滲流連續(xù)性方程推導(dǎo)的v_z，右側(cè)代表基于進氣值面處水量均衡關(guān)系推導(dǎo)的v_z），共8個進氣值面處v_z的解析表達進行討論。

3.2.1"線性Boussinesq方程進氣值面處的v_z

1）完整解計算的進氣值面處v_z

仿Theis井流模型采用了線性的定位在進氣值面處的Boussinesq方程（式（5a）），該方程等號右側(cè)是基于進氣值面處水量均衡關(guān)系推導(dǎo)的v_z，記作v_z，ha1，等號左側(cè)是基于進氣值面處滲流連續(xù)性方程推導(dǎo)的v_z，記作v_z，ha2。當(dāng)采用極坐標(biāo)表示，并將式（5a）中的s′_ha用式（4）替換成h時，則v_z，ha1和v_z，ha2如下式所示：

v_z，ha1=μ_haht； """"""""""""（21a）v_z，ha2=Kh_cp²hr²+Kh_cprhr+Kh_cph（hr）²"。（21b）

將進氣值面高程的完整解式（9）代入式（21），則v_z，ha1和v_z，ha2可分別寫為：

v_z，ha1=－Qμ_hae^－r2μ_ha4Tt4πKth₀²－Q∫_r2μha4Tt1ye^－ydy2πK¹²"；（22a）v_z，ha2=－Qμ_hae^－r2μ_ha4Tt4πKth₀²－Q∫_r2μha4Tt1ye^－ydy2πK¹²。（22b）

由式（22）可知，v_z，ha1和v_z，ha2結(jié)果非0且相等，表明基于定位在進氣值面處的Boussinesq方程可以計算得到v_z。

將表3獲得的參數(shù)反演結(jié)果代入式（22）中，計算得到a、b、c井抽水試驗對應(yīng)的進氣值面處v_z，計算結(jié)果如圖4所示。由圖4可見，當(dāng)三口不同位置（a、b、c井，r=2.60、1.90、3.50 m）的抽水井定流量抽水時，隨抽水時間的延長，進氣值面處的v_z呈非線性增長的趨勢。在抽水試驗早期，滲流分速度急速上升；隨抽水時長增加，滲流分速度上升速率減小，進氣值面處的v_z趨于平緩。

2）近似解計算的進氣值面處v_z

將進氣值面高程的近似解式（12）代入式（21），則基于進氣值面處水量均衡關(guān)系推導(dǎo)的v_z記作v′_z，ha1，基于進氣值面處滲流連續(xù)性方程推導(dǎo)的v_z記作v′_z，ha2，可分別寫為：

v′_z，ha1=－Qμ_ha4πKt（h₀²－

Q2πK

ln2.25Ttr²μ_ha）^－12;（23a）v′_z，ha2=0 。""""""""""""（23b）

由式（23）可知，v′_z，ha1結(jié)果非0，v′_z，ha2結(jié)果為0，因此可以用基于進氣值面處水量均衡關(guān)系來計算v_z。將表3獲得的參數(shù)反演結(jié)果代入式（23a）中計算得到a、b、c井抽水試驗對應(yīng)的進氣值面處的v_z，計算結(jié)果如圖4所示。

3.2.2"非線性Boussinesq方程進氣值面處的v_z

1）完整解計算的進氣值面處v_z

對于非線性處理的定位在進氣值面處的Boussinesq方程，等號右側(cè)是基于進氣值面處水量均衡關(guān)系推導(dǎo)的v_z，記作v″_z，ha1，等號左側(cè)是基于進氣值面處滲流連續(xù)性方程推導(dǎo)的v_z，記作v″_z，ha2：

v″_z，ha1=μ_haht ；""""""""""（24a）v″_z，ha2=Kx（hhx）+y（hhy） 。"（24b）

將式（24）轉(zhuǎn)化成極坐標(biāo)下的表達式之后，再將進氣值面高程的完整解（式（9））代入式（24），即獲得基于非線性方程的v_z近似表達：

v″_z，ha1=－Qμ_hae^－r2μ_ha4Tt4πKth₀²－Q∫_r2μha4Tt1ye^－ydy2πK¹²；"""""""""""""""""""（25a）v″_z，ha2=－Q²e^－r2μ_ha2Tt4π²Kr²h₀²－Q∫_r2μha4Tt1ye^－ydy2πK¹²－""""Qμ_hae^－r2μ_ha4Tt4πTt。""""""""（25b）

由式（25）可知，利用進氣值面高程的完整解代入非線性處理的定位在進氣值面處的Boussinesq方程的等號兩側(cè)計算得到的v″_z，ha1和v″_z，ha2非0且不相等。

將表3獲得的參數(shù)反演結(jié)果代入式（25）中計算得到a、b、c井抽水試驗對應(yīng)的進氣值面處v_z，計算結(jié)果如圖4所示。由圖4可見，當(dāng)三口不同位置（a、b、c井，r=2.60、1.90、3.50 m）的抽水井定流量抽水時，完整解計算的進氣值面處v_z的兩條曲線均表現(xiàn)為隨時間的增長呈非線性增長的趨勢。在抽水試驗早期，滲流分速度急劇增長；隨抽水時長增加，滲流分速度上升速率減小，進氣值面處的v_z趨于平緩。兩者的變化趨勢是相似的，但也存在一定的差異。

2）近似解計算的進氣值面處v_z

將進氣值面高程的近似解（式（12））代入式（24），還可獲得基于進氣值面處水量均衡關(guān)系推導(dǎo)v_z近似解的另一種表述，記作v_z，ha1，以及基于進氣值面處滲流連續(xù)性方程推導(dǎo)v_z的另一種表述，記作v_z，ha2，可分別寫為：

v_z，ha1=－Qμ_ha4πKt（h₀²－

Q2πK

ln2.25Ttr²μ_ha）^－12";（26a）v_z，ha2=－Q²4π²Kr²（h₀²－Q2πKln2.25Ttr²μ_ha）－12。

（26b）

由式（26）可知，利用進氣值面高程近似解代入非線性化處理的定位在進氣值面處Boussinesq方程的等號兩側(cè)計算得到的v_z，ha1和v_z，ha2非0且不相等，兩者之間相差μ_har²π/（Qt）倍。

將表3獲得的參數(shù)反演結(jié)果代入式（26）中計算得到a、b、c井抽水試驗對應(yīng)的進氣值面處v_z，計算結(jié)果如圖4所示。由圖4可見，當(dāng)三口不同位置（a、b、c井，r=2.60、1.90、3.50 m）的抽水井定流量抽水時，v_z，ha1隨抽水時間的延長呈非線性增長的趨勢，且在抽水試驗早期，滲流分速度急速上升；隨抽水時長增加，滲流分速度上升速率減小，進氣值面處的v_z趨于平緩。v_z，ha2隨抽水時間的延長呈下降趨勢，但下降速率不大，基本趨于平緩。

3.2.3"不同方法計算的進氣值面處v_z的對比

由圖4可知，式（22a）（v_z，ha1）、式（22b）（v_z，ha2）、式（25a）（v″_z，ha1），式（25b）（v″_z，ha2），式（23a）（v′_z，ha1）、式（26a）（v_z，ha1）計算得到的進氣值面處v_z的趨勢相似，均隨時間的增長呈非線性增長趨勢。其中，v_z，ha1、v_z，ha2、v″_z，ha1相等，v′_z，ha1、v_z，ha1與v_z，ha1、v_z，ha2、v″_z，ha1基本相等，而v″_z，ha2與v_z，ha1、v_z，ha2、v″_z，ha1差異較大。導(dǎo)致v″_z，ha2出現(xiàn)較大偏差的原因在于，將基于線性Boussinesq方程計算獲得的進氣值面高程近似為非線Boussinesq方程計算的進氣值面高程，計算得到的v_z會出現(xiàn)較大偏差。若在式（25b）的基礎(chǔ)上，將利用式（10）近似函數(shù)W（u_μ）獲取的進氣值面高程近似解替換進氣值高程完整解，即式（25b）變換為式（26b），則該式計算獲得進氣值面處v_z隨時間發(fā)展的變化趨勢（v_z，ha2）與式（25b）計算獲得變化趨勢（v″_z，ha2）相反；表現(xiàn)為進氣值面處v_z隨時間的增長呈非線性遞減趨勢，該趨勢也與式（22a）、式（22b）、式（25a）的計算獲得的變化趨勢相反。因此，不推薦使用式（25b）和式（26b）計算進氣值面處v_z。

由式（23b）計算得到的進氣值面處v_z恒為0。導(dǎo)致該結(jié)果的原因是利用式（10）對W（u_μ）函數(shù)近似。該結(jié)果不符合完整井定流量抽水時存在鉛直流分速度的基本認(rèn)識，故而不推薦使用式（23b）計算進氣值面處的v_z。

綜上可知，無論進氣值面處的Boussinesq方程是否經(jīng)過線性處理，利用完整解（式（9））和近似解（式（12））分別計算得到的進氣值面處v_z（式（22a）、式（22b）、式（25a）、式（23a）、式（26a））的數(shù)值基本吻合，曲線趨勢也一致。建議在實際應(yīng)用中，可采用上述公式計算進氣值面處的v_z。不推薦使用以下三個公式：近似解代入線性Boussinesq方程中的基于進氣值面處滲流連續(xù)性方程推導(dǎo)的v_z（式（23b））、完整解或近似解代入非線性Boussinesq方程中基于進氣值面處滲流連續(xù)性方程推導(dǎo)的v_z（式（25b）和式（26b））。

3.3"優(yōu)點與不足

本文采用定位在進氣值面處的Boussinesq方程作為控制方程修正仿Theis井流模型，并獲取修正模型的解析解，豐富和發(fā)展了非穩(wěn)定井流相關(guān)理論。本文提出的修正仿Theis井流模型是對經(jīng)典仿Theis井流模型的修正。根據(jù)仿Theis井流模型的具體應(yīng)用可以預(yù)測本文的修正仿Theis井流模型的應(yīng)用前景。仿Theis井流公式被廣泛應(yīng)用到水文地質(zhì)、工程地質(zhì)中，如估計水文地質(zhì)參數(shù)，預(yù)測水位降深、出水量和影響半徑等。修正仿Theis井流除了與經(jīng)典仿Theis井流模型具有相似的應(yīng)用場景之外，還可以與熱傳輸方程耦合，模擬潛水抽水的水熱耦合傳輸過程，或與溶質(zhì)運移方程耦合，模擬潛水抽水驅(qū)動下的溶質(zhì)運移過程，等等。然而，本文方法仍存在以下不足：

1）本文的理論基礎(chǔ)之一是土水特征曲線的Brooks-Corey模型^[29]。該模型忽略了滲流演化歷史對土水特征曲線的影響，將土水特征曲線簡化為單值性曲線（即含水率與毛細(xì)壓強水頭一一對應(yīng)），對于給定的均質(zhì)沉積物其進氣值唯一確定。然而，進氣值面高度變化是復(fù)雜的物理過程，其受到地層介質(zhì)的非均質(zhì)性、地層結(jié)構(gòu)特征、地下水下降速度等多方面因素的影響。在后繼的研究中應(yīng)結(jié)合上述因素影響下進氣值特性的研究成果，進一步完善修正仿Theis井流模型，拓展模型的適用范圍，提高模型對滲流過程刻畫的精度。

2）本文提出的修正仿Theis井流模型，建立在定位在進氣值面處的Dupuit假定基礎(chǔ)之上，與經(jīng)典仿Theis井流模型（建立在定位在潛水面處的Dupuit假定基礎(chǔ)之上）具有相似的局限性，即只適用于坡度較平緩、等水頭面近似鉛直的情形，適合計算小降深或抽水后期的降深。

4"結(jié)論

1）修正經(jīng)典仿Theis井流模型，將自由面由潛水面移動到進氣值面，采用定位在進氣值面處的Boussinesq方程作為修正的仿Theis井流模型控制方程，通過線性化方法和拉普拉斯變換，推導(dǎo)出修正井流模型的解析解，并基于新的解析解構(gòu)建參數(shù)反演模型，拓展了仿Theis模型。

2）通過將修正模型計算的降深與實測降深進行對比，結(jié)果表明修正的仿Theis井流模型計算的潛水面處降深與實測降深基本一致，驗證了本文所提出的井流模型合理性。

3）對修正的仿Theis井流模型進行參數(shù)敏感性分析表明，進氣值面高程隨與井中心的距離和給水度增大而增大，隨抽水流量增大而減小。在抽水早期，進氣值面高程隨滲透系數(shù)的增大而減小，在抽水后期則相反。

4）將進氣值面高程的完整解和近似解代入（非）線性處理的定位在進氣值面處的Boussinesq方程，共獲取8個進氣值面處鉛直方向滲流分速度的解析表達。其中，對于線性處理的定位在進氣值面處的Boussinesq方程，將進氣值面高程的完整解代入基于進氣值面處水量均衡關(guān)系和基于進氣值面處滲流連續(xù)性方程推導(dǎo)的鉛直方向滲流分速度的表達式，結(jié)果非0且相同。利用這8個解析表達計算完整井定流量抽水條件下的進氣值面處鉛直方向滲流分速度的計算結(jié)果表明：基于進氣值面處水量均衡關(guān)系（基于線性或非線性方程，采用完整解或近似解）以及基于進氣值面處滲流連續(xù)性方程（基于線性方程、采用完整解）等5種情況下推導(dǎo)的鉛直方向滲流分速度的解析表達計算結(jié)果基本吻合，曲線趨勢一致，表現(xiàn)為隨抽水時間的延長呈非線性增長。

5）除此之外，對于線性處理的定位在進氣值面處的Boussinesq方程，將完整解代入基于進氣值面處滲流連續(xù)性方程推導(dǎo)的鉛直方向滲流分速度的結(jié)果恒為0；對于非線性處理的定位在進氣值面處的Boussinesq方程，將完整解或近似解代入基于進氣值面處滲流連續(xù)性方程推導(dǎo)的鉛直方向滲流分速度與其余5種情況均存在較大誤差。

參考文獻（References）：

[1] 王玉林，謝康和，李傳勛，等. 抽-灌同軸非完整井承壓層非穩(wěn)定流模型及解析解[J].水利學(xué)報，2012，43（1）：60-68.

Wang Yulin， Xie Kanghe， Li Chuanxun， et al. A Mathematical Model and Its Analytical Solution for Cconfined Aquifer Subjected to Pumping and Recharge Implemented by a Single Partially Penetrating Well[J]. Journal of Hydraulic Engineering， 2012， 43（1）：60-68.

[2] 李文鵬.對地表水資源與地下水資源重復(fù)量的認(rèn)識與水資源開發(fā)利用理念探討[J]. 水文地質(zhì)工程地質(zhì)，2023，50（1）：1-2.

Li Wenpeng. Discussion on the Recognition of Surface Water Resources and Groundwater Resources Duplication and the Concept of Water Resources Development and Utilization[J]. Hydrogeology amp; Engineering Geology， 2023，50（1）：1-2.

[3] 趙寶峰， 康衛(wèi)東， 馬蓮凈，等. 抽水試驗和長觀水位聯(lián)合模擬確定含水層參數(shù)[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版），2009，39（3）：482-486.

Zhao Baofeng， Kang Weidong， Ma Lianjing，et al. Aquifer Parameter Recognition by Combining Simulation of Pumping Test and Water Level of Long-Term Observation Well[J]. Journal of Jilin University （Earth Science Edition）， 2009，39（3）：482-486.

[4] 陳晨. 承壓含水層中變流量抽水試驗井流動力學(xué)模型與參數(shù)反演方法研究[D].武漢：中國地質(zhì)大學(xué)，2022. Chen Chen. Study on Flow Dynamics Model and Parameter Inversion Method of Variable Flow Pumping Test Well in Confined Aquifer[D]. Wuhan： China University of Geosciences，2022.

[5] 李丞，王寧濤，胡成，等. 利用多孔抽水試驗確定承壓-無壓流條件下承壓含水層的水文地質(zhì)參數(shù)[J].安全與環(huán)境工程，2020，27（4）：1-7.

Li Cheng， Wang Ningtao， Hu Cheng， et al. Hydrogeological Parameter Determination of the Confined Aquifer Under Confined-Unconfined Conditions by Means of Multiple-Well Pumping Test[J]. Safety and Environmental Engineering， 2020，27（4）：1-7.

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

Qiu Shuwei， Wu Yamin， Ke Yuqi， et al. Optimization of Hydrogeological Parameters Based on Ergodic Search Algorithm[J]. Journal of Jilin University （Earth Science Edition）， 2020，50（6）：1854-1861.

[7] Theis C V. The Relation Between the Lowering of the Piezometric Surface and the Rate and Duration of Discharge of a Well Using Ground Water Storage[J]. American Geophysical Unions Transactions， 1935，16（2）： 519-524.

[8] Yeh H D， Yang S Y， Peng H Y. A New Closed-Form Solution for a Radial Two-Layer Drawdown Equation for Groundwater Under Constant-Flux Pumping in a Finite-Radius Well[J]. Advances in Water Resources ，2003，26（7）：747-757.

[9] Yeh H D， Chang Y C. Recent Advances in Modeling of Well Hydraulics[J]. Advances in Water Resources， 2013，51：27-51.

[10] Wen Z， Zhan H B， Wang Q R， et al. Well Hydraulics in Pumping Tests with Exponentially Decayed Rates of Abstraction in Confined Aquifers[J]. Journal of Hydrology， 2017，548：40-45.

[11] Smith R， Li J. Modeling Elastic and Inelastic Pumping-Induced Deformation with Incomplete Water Leveal Records in Parowan Valley，Utah[J]. Journal of Hydrology，2021，601：126654.

[12] 袁帥，王君，吳朝峰，等.虹吸排水法處理軟土地基的水位與沉降計算模型[J].吉林大學(xué)學(xué)報（地球科學(xué)版）， 2024，54（1）：208-218.

Yuan Shuai，Wang Jun，Wu Zhaofeng，et al. Calculation Model for Water Level and Settlement of Soft Foundation Treated by Siphon Drainage[J]. Journal of Jilin University （Earth Science Edition），2024， 54（1）：208-218.

[13] Mehrem M F A A. Fluid-Dynamic and Geomechanical Modelling of the Shallow Formations of the PO Plain[D]. Turin： Politecnico di Torino，2022.

[14] 郭東屏，宋焱勛，錢會，等.地下水動力學(xué)[M].西安：陜西科學(xué)技術(shù)出版社， 1994.

Guo Dongping，Song Yanxun，Qian Hui，et al. Hydraulics of Groundwater[M]. Xi’an： Shaanxi Science and Technology Press，1994.

[15] 鄧厚基.試算法借助輔助曲線圖一次確定u₁值簡化計算a、T（K）的方法[J].河北地質(zhì)學(xué)院學(xué)報，1978（3）：80-83.

Deng Houji. The Experimental Algorithm Simplifies the Method of Calculating a and T（K） by Determining the u₁"Value at Once with the Aid of the Auxiliary Graph[J]. Journal of Hebei College of Geology，1978（3）：80-83.

[16] 孫德全，魯孟勝，張兆民.新疆大南湖北露天煤礦首采區(qū)Ⅲ火燒區(qū)地下水資源的數(shù)值模擬[J].煤田地質(zhì)與勘探，2014，42（4）：64-68.

Sun Dequan， Lu Mengsheng， Zhang Zhaomin.The Numerical Simulation of Groundwater Resources in Burnt Zone of the First Mining Area Ⅲ in Dananhu Northern Surface Mine of Xinjiang[J]. Coal Geology amp; Exploration，2014，42（4）：64-68.

[17] 胡冰.石家莊無極縣某化工廠項目地下水環(huán)境影響及預(yù)測研究[D].長春：吉林大學(xué)，2018.

Hu Bing. Prediction and Evalution of Groundwater Contamination for a Chemical Plant Project in Wuji County， Shijiazhuang[D]. Changchun： Jilin University，2018.

[18] 王曉惠，束龍倉，鄧銘江，等.非穩(wěn)定流非完整輻射井抽水試驗水文地質(zhì)參數(shù)的求解[J].水電能源科學(xué)，2011，29（7）：43-46，114.

Wang Xiaohui， Shu Longcang， Deng Mingjiang， et al.Hydrogeologic Parameters of Pumping Test for Incomplete Radial Well Under Unsteady Flow[J]. Water Resources and Power， 2011，29（7）：43-46，114.

[19] 王俊杰，張慧萍.地下水水位預(yù)測在深基坑工程中的應(yīng)用[J].地質(zhì)災(zāi)害與環(huán)境保護，1999，10（4）：63-66.

Wang Junjie， Zhang Huiping. The Application of Groundwater Table Predication in Deep Foundation Excavation[J]. Journal of Geological Hazards and Environment Preservation， 1999，10（4）：63-66.

[20] 許家美，谷鴻飛，俞裕香.泰斯公式在軟土地區(qū)管井降水中的應(yīng)用[J].地下水，2007，29（1）：121-124.

Xu Jiamei， Gu Hongfei， Yu Yuxiang. Theis Formula’s Application of Lowering Water by Tube Well in the Soft-Land Area[J]. Ground Water， 2007，29（1）：121-124.

[21] Boulton N S. The Influence of Delayed Drainage on Datafrom Pumping Tests in Unconfined Aquifers[J]. Journal of Hydrology， 1973，19（2）：157-169.

[22] Neuman S P. Theory of Flow in Unconfined Aquifers Considering Delayed Response of the Water Table[J]. Water Resources Research， 1972，8（4）：1031-1045.

[23] 陳崇希，唐仲華.Theis不穩(wěn)定潛水井流模型的改進：具入滲補給[J].水文地質(zhì)工程地質(zhì)，2021，48（6）：1-12.

Chen Chongxi， Tang Zhonghua. Improvement of the Theis Unsteady Well Flow Model with Infiltration Recharge in a Phreatic Aquifer[J]. Hydrogeology amp; Engineering Geology， 2021，48（6）：1-12.

[24] Bear J. Hydraulics of Groundwater[M]. McGraw-Hill Series in Water Resources and Environmental Engineering[M]. New York：McGraw-Hill International Book Company， 1979.

[25] Cheng D W，Zhan H B，Li J ， et al. Characteristics of Inverted Saturated Zone Under Unclogged Streams[J]. Journal of Hydrology， 2021，597（4）：126288.

[26] 程大偉.定位在進氣值面的Boussinesq方程[C]//第八屆青年地學(xué)論壇.武漢：長江科學(xué)院院報，2023.

Cheng Dawei. The Boussinesq Equation Located at the Air Entry Plane[C]//The 8th Young Scientist Forum of Earth Science. Wuhan： Journal of Changjiang River Scientific Research Institute， 2023.

[27] 江峰，黃琨，馬丁山，等.潛水含水層井流模型的對比研究[J].地質(zhì)科技情報，2015，34（3）：160-164.

Jiang Feng， Huang Kun， Ma Dingshan， et al.Comparative Study on Well Flow Model in Unconfined Aquifer[J]. Geological Science and Technology Information， 2015，34（3）：160-164.

[28] 朱悅璐，熊俊飛，鄧煒.基于分段函數(shù)的潛水井降深曲線在工程實踐中的應(yīng)用[J].人民黃河，2022，44（6）：70-75.

Zhu Yuelu， Xiong Junfei， Deng Wei.Application of Phreatic Well Drawdown Curve Based on Piecewise Function in Engineering[J]. Yellow River， 2022，44（6）：70-75.

[29] Brooks R H， Corey A T. Hydraulic Properties of Porous Media and Their Relation to Drainage Design[J]. Transactions of the ASAE， 1964， 7（1）：26-28.