方神光
(珠江水利委員會 珠江水利科學(xué)研究院, 廣東 廣州 510611)
珠江河口屬典型弱潮型河口, 一定條件下, 容易出現(xiàn)海水與河水之間明顯的分層現(xiàn)象[1], 海水在表層淡水下方呈楔狀沿河底向上游推進, 形成鹽水楔或異重流。如在伶仃洋水域, 來自南海大陸架高鹽水體在潮汐作用力下, 主要沿伶仃水道南段和暗士頓水道上溯[2], 與虎門口下泄徑流形成混合和分層現(xiàn)象。珠江河口八大口門中, 以磨刀門輸水和輸沙量最大, 磨刀門水道上游是中山、珠海和澳門的主要水源地, 近年由于氣候變化及人類活動等, 河口咸潮入侵頻繁, 嚴重影響沿岸城市的供水安全。
鑒于磨刀門水道的重要性, 有關(guān)該水道咸潮入侵的原因和機理研究成果較多, 如劉杰斌等[3]和包蕓等[4]分析實測資料顯示, 磨刀門水道出現(xiàn)枯季小潮期間咸潮快速上溯, 是由小潮期間磨刀門水道凈泄量為負所致; 聞平等[5]認為影響磨刀門水道咸潮入侵的主要因素是徑流, 并推薦了最小壓咸流量范圍以及最佳壓咸補淡的時機; 陳榮力等[6]和羅丹[7]分析實測資料, 認為大潮期間加大上游下泄流量是最好的壓咸時機; 韓志遠等[8]的研究顯示, 近些年的河道采砂以及河口圍墾導(dǎo)致磨刀門口門“調(diào)淡”作用的喪失, 才是咸潮影響加劇的主要原因; 河口咸潮入侵與徑流和潮汐作用的此消彼長關(guān)系密切, 賈良文等[9]的分析顯示, 磨刀門水道枯季潮流為主要動力,河口下層有反向水流, 存在明顯的因密度差而形成的密度環(huán)流。在磨刀門口, 當(dāng)上游西江梧州來流量小于2 500 m3/s時, 屬于以潮為主河口, 潮流動力作用加強, 咸潮對上游的影響逐步加大。以西蒙斯咸淡水混合判別指標(biāo)分析, 咸淡水混合狀態(tài)多為弱混合型及緩混合型, 屬高度分層型, 形成明顯鹽水楔[10]。
可見, 枯季磨刀門水道由于徑流減少, 鹽水入侵容易形成分層現(xiàn)象, 具有典型的三維水動力特性。鑒于此, 依托磨刀門水道最新地形資料, 本文采用Delft3D軟件建立了磨刀門水道三維潮流和鹽度數(shù)值模型, 并采用2009年12月10~25日的實測資料對模型進行了率定和驗證。在此基礎(chǔ)上, 對磨刀門水道的潮流和鹽度入侵特性進行了分析和探討。
控制方程由連續(xù)方程、動量方程和鹽度方程組成, 在 Delft3D中, 采用曲線貼體σ坐標(biāo)系, 所建立的數(shù)值模型中已考慮鹽度正壓和斜壓對水流運動的影響, 以便更真實地模擬河口徑流和潮汐的相互作用。其通用形式如公式(1)。
其中,φ表示流場中不同物理量的通用變量,Rφ(xi)表示源項。uj是貼體坐標(biāo)系上xj的速度分量;ρ是水體密度, 隨溫度和鹽度變化, 采用國際通用海水密度計算公式(UNESCO), 模型中已考慮由鹽度變化引起的水體密度梯度對水流運動的影響;Γφ為擴散項參數(shù)。DELFT3D軟件中的相關(guān)控制方程詳細表達形式可參考相關(guān)文獻[11-12], 并采用交替方法對控制方程進行離散和求解。
圖1 磨刀門水道實測站點分布圖Fig.1 Sketch of field station locations in Modaomen Waterway
2009年 12月10~25日(農(nóng)歷 10月 24日~11月10日), 珠江水利科學(xué)研究院[7]對磨刀門水道竹銀-磨刀門口門(大橫琴水文站以南約3.9 km)之間近40 km的磨刀門水道中的 8個站點的潮位、流速和鹽度進行了逐時測量, 實測點位置如圖1所示, 其中1#測點位于磨刀門口門位置, 距大橫琴水文站以南約3.8 km,8#站點位于竹銀, 4#站點布置在洪灣水道內(nèi)。數(shù)值模擬試驗中, 為方便驗證和提取邊界, 選取的模型計算區(qū)域與實測區(qū)域一致(圖1)。采用正交曲線網(wǎng)格對計算區(qū)域進行剖分, 共計 216×27個網(wǎng)格, 垂直方向平均分成10層。
1.3.1 邊界條件
計算范圍內(nèi)的磨刀門水道邊界條件相對較為簡單, 以該水道上游竹銀位置為上邊界, 應(yīng)用2009年12月10~25日實測流量數(shù)據(jù)進行控制, 下邊界主要是磨刀門河口邊界和洪灣水道邊界, 采用同一時期實測潮位控制。
1.3.2 初始條件
DELFT3D 軟件中, 初始條件有“冷啟動”和“熱啟動”之分, “冷啟動”即給定一個零初始場, 依靠模型不斷迭代計算以接近實際情況; “熱啟動”即給定初始水流、水位、鹽度等初始場。在潮位和潮流的數(shù)值模擬計算中, 邊界影響能迅速傳遞到流場內(nèi)部,因此模擬的流態(tài)能迅速趨向于真實值; 但在鹽度場計算中, 由于鹽度輸移和擴散緩慢以及遵循物質(zhì)守恒規(guī)律, 鹽度若給為零初始場, 則要花費相當(dāng)長的計算周期才能達到真實分布狀態(tài)。因此, 如果給定的鹽度場能盡量接近真實情況, 能提高收斂速度, 節(jié)省計算時間。鑒于此, 本文依據(jù)2009年枯季磨刀門水道布置的 8個實測點位的鹽度資料以及實際地形狀況, 采用插值方法給出了計算區(qū)域每一層網(wǎng)格節(jié)點上的鹽度初始值。
1.3.3 其他參數(shù)選取
根據(jù)模型調(diào)試顯示, 模型糙率變化范圍為0.016~0.025, 從上游至口門依次減少。由于水平方向相對垂直方向大的多, 因此水平方向的紊動黏性系數(shù)和鹽度擴散系數(shù)根據(jù)模型調(diào)試后的結(jié)果, 直接給定為常值即可; 在垂直方向, 則通過水動力k-ε(紊動能-紊動耗散率)紊流模型計算給出。
由于選取1#、8#和4#站點分別作為磨刀門口門、磨刀門水道上游以及洪灣水道的邊界控制站點, 因此該3個測站不參與驗證。圖2~圖6給出了2#和6#測站的潮位、表底層流速和鹽度的驗證圖, 其他測站由于篇幅所限不一一列出。驗證結(jié)果顯示, 從潮位、流速和鹽度各項實測值與數(shù)值模擬結(jié)果比較來看,數(shù)值模擬計算的潮位和流速與實測值吻合較好, 潮位的絕對誤差在 0.07 m以內(nèi), 其中 2#的誤差最小,6#站的誤差最大。流速除了個別站位存在較大的誤差, 誤差均在10%左右。鹽度吻合相對較差, 個別實測點底層鹽度數(shù)值模擬結(jié)果與實測誤差較大, 尤其是底層靠近上游的驗證站位, 這可能是由該測量站位資料存在一定誤差引起。但從總體來看, 模擬計算各測點鹽度變化總體趨勢與實測基本一致, 綜合考慮和分析實測過程中各種干擾因素對鹽度測量值的影響、局部地形的差異以及數(shù)值模擬的誤差精度等,本模型基本能反映 2009年 12月 10~25日(農(nóng)歷 10月24日~11月10日)期間的潮位、流速和鹽度的變化, 可以應(yīng)用本模型來分析該水道潮流和鹽度變化規(guī)律。
圖2 潮位驗證圖Fig.2 Verification of tidal level
圖3 2#站流速驗證圖Fig.3 Verification of velocities in No.2 station
圖4 6#站流速驗證圖Fig.4 Verification of velocities in No.6 station
圖5 2#站鹽度驗證圖Fig.5 Verification of salinities in No.2 station
圖6 6#站鹽度驗證圖Fig.6 Verification of salinities in No.6 station
圖7給出了數(shù)值模擬的大潮漲落急時刻表層和底層的平面流態(tài)(中層及中潮和小潮時刻流態(tài)與大潮基本一致, 此處沒有給出), 圖8和圖9給出了實測情況下的表層和底層流速沿河口向上游方向的變化線,橫坐標(biāo)表示以口門 1#測點為起始原點, 向上游方向距離為正。從平面流態(tài)來看, 靠近東側(cè)為磨刀門水道主槽, 西側(cè)主要為淺灘分布, 因此東側(cè)流速要明顯大于西側(cè)。大潮漲急時刻, 從口外進入磨刀門水道的漲潮流向北, 即西北方向; 洪灣水道水流為東南向,即磨刀門水道部分漲潮流進入到該水道。從表、底分層來看, 表層平均流速略大于底層平均流速。大潮落急時刻, 落潮流速顯著大于漲急流速, 磨刀門水道落潮流為東南向; 洪灣水道流向則為西北向, 部分潮流進入磨刀門水道。流速存在分層現(xiàn)象, 由表向底流速依次減小。對中潮和小潮的分析顯示, 磨刀門水道表層落潮流速顯著大于漲潮流速, 大潮和中潮時, 底層落潮流速略大于漲潮流速, 小潮時, 底層落潮流速略小于漲潮流速; 磨刀門水道表層平均漲潮流速值隨潮型增大而增大, 顯示漲潮時潮流在磨刀門水道內(nèi)占優(yōu)勢; 落潮流時, 徑流作用凸顯, 與落潮流相互作用, 導(dǎo)致河道流速變化較為復(fù)雜, 隨潮型變化趨勢不明顯; 各潮型下, 底層漲、落潮流速都隨潮型增大而增大。總體來看, 枯季磨刀門水道總體漲、落潮流速都不大, 各潮型下的表層總體漲潮平均流速都在0.5 m/s以內(nèi), 總體落潮平均流速在0.8 m/s以內(nèi); 底層總體漲落潮平均流速都在0.5 m/s以內(nèi)。
枯季小潮期間, 徑流量小, 口門潮汐動力弱, 容易出現(xiàn)鹽水入侵造成的密度分層現(xiàn)象, 表層淡水向海漂移, 而底層鹽水由于補償流和鹽度斜壓等作用,向上游運動。圖10給出了小潮時刻(12月 11日 01時小潮低潮)磨刀門水道的縱剖面流速等值線, 此時流速處于漲落潮的交替時刻??v坐標(biāo)為高程, 坐標(biāo)原點設(shè)為口門位置 1#測點, 橫坐標(biāo)為以口門 1#測點為原點向上游方向的距離。圖中流速為正表示流向向海, 負值表示向陸地??梢? 小潮低潮時刻, 底層水流以向陸為主, 表層水流以向海運動為主, 表底層流向相反; 垂向上, 流向變化的位置流速等值線分布較為密集, 流速梯度變化大。實測期間, 該現(xiàn)象在口門 1#測點能明顯觀察到。大、中、小潮型下, 1#測點表層平均落潮和漲潮歷時分別為8.8 h和3.9 h、12.5 h和3.8 h、8.6 h和3.9 h; 底層平均落潮和漲潮歷時分別為8 h和16.8 h、7.3 h和17 h、6 h和18.9 h。該測點表層和底層漲落潮時段不完全重合, 說明表層和底層流向存在不一致的時段; 從各潮型來看,表、底層流速相反的時間段一般出現(xiàn)在經(jīng)過第一個長時間的落潮流之后的第二個較短落潮時間段內(nèi)。且沿水深方向流向發(fā)生變化的位置隨潮動力減弱而有向水面移動的趨勢。
圖7 磨刀門水道大潮期間漲落急流態(tài)Fig.7 Flowing at Modaomen waterway during spring
圖10 磨刀門水道小潮縱剖面流速圖Fig.10 Velocity contour at a longitudinal cross section during neap
圖11給出了數(shù)值模式計算得到的大潮漲落急時刻表層和底層的鹽度分布, 圖12給出了觀測站點表層和底層鹽度分布。由于磨刀門水道從口門至與洪灣水道匯合口之間河段的主槽靠近東岸, 因此表底層鹽度入侵受地形影響明顯, 總體呈現(xiàn): 漲潮時磨刀門水道東側(cè)鹽度大于西側(cè), 落潮時東側(cè)鹽度小于西側(cè)的趨勢。漲急時, 表層和底層鹽度為10的等值線分別到達掛定角和大沖口水閘附近; 落急時, 表層和底層鹽度為10的等值線線則都到達燈籠水閘附近??梢? 漲潮時, 鹽度為 10的等值線入侵距離超過表層; 落潮時, 兩者相差不大; 實測和模擬結(jié)果都顯示, 磨刀門水道表層鹽度顯著小于底層鹽度,小潮期間該現(xiàn)象最為明顯, 且小潮期間的底層鹽度要大于大潮和中潮, 顯示磨刀門水道枯季小潮期間的鹽水入侵更為嚴重, 并形成明顯的鹽水分層現(xiàn)象,而表層鹽度變化規(guī)律仍基本遵循隨潮型增大而增大的規(guī)律。竹排沙淺灘東側(cè)水道鹽度總體小于西側(cè), 大排沙淺灘東側(cè)石岐水道鹽度則明顯大于西側(cè)磨刀門水道。從表、底層鹽度等值線分布來看, 漲落潮時底層鹽度等值線分布梯度顯著大于表層, 同一位置的底層鹽度明顯高于表層。
為探討大、中、小潮下磨刀門水道鹽度在垂向的變化規(guī)律, 圖13給出沿主槽對應(yīng)剖面漲落潮下的鹽度垂向分布, 縱剖面線位置如圖1所示。由于磨刀門水道地形變化對鹽水入侵影響明顯, 為分析方便,根據(jù)水深地形特征, 大致可以將磨刀門水道近口門河段分為4段: 河段1為從口門至與洪灣水道匯合口,河段長度大致為15 km, 從口門向河道上游方向14 km范圍內(nèi), 河底地形高程從大約–4 m下降至接近–12 m,然后在1 km內(nèi)迅速上升到約–3 m; 河段2為從洪灣水道匯合口到竹排沙, 該河段長約16 km, 呈現(xiàn)兩端高、中間低的形態(tài), 該河段下游端地形高程為–3 m,上游端地形高程接近 0, 中間最低位置地形高程為–12 m左右; 河段3為石岐水道, 該水道位于大排沙的東側(cè), 河道長度約6 km, 從縱剖面地形上看, 水深較淺, 也呈現(xiàn)兩端高、中間低的形式; 河段4為從磨刀門水道進入石岐水道分流口往上游方向河段, 該河段河床地形由–2 m高程在5 km距離內(nèi)迅速下降到–15 m。根據(jù)各潮型下的縱剖面鹽度分析各河段鹽度特性可見:
圖11 磨刀門水道大潮期間鹽度分布Fig.11 Salinity contours at Modaomen waterway during spring
1) 大潮漲急時, 外海高鹽度(鹽度在20以上)的純海水進入到磨刀門水道河段 1內(nèi), 表層和底層入侵距離分別達4 km和10 km左右, 底層高鹽水體入侵距離顯著大于表層, 高鹽海水與沖淡水接觸鋒面鹽度等值線分布密集, 鹽度等值線近似垂向分布,與洪灣水道匯合口處鹽度約為 12; 河段 2底層鹽水入侵距離大于表層, 末端竹排沙位置鹽度為 4左右;河段3末端西河水閘位置鹽度約為1; 河段4從西河水閘向上游, 鹽度都小于1。
2) 大潮落急時, 河段 1表層高鹽海水基本退出口門, 為沖淡水占據(jù), 底層由于地形原因, 部分高鹽水體聚集在底部無法退出, 形成顯著的表底分層現(xiàn)象, 與洪灣水道匯合口處鹽度為 16左右, 明顯高于漲潮時的鹽度; 河段2表、底層鹽度向上游入侵的距離反而大于大潮漲急時, 該河段末端竹排沙位置鹽度為 6左右; 河段3整體鹽度較漲急時有顯著升高,末端西河水閘位置鹽度約為5; 河段4的整體鹽度也呈現(xiàn)顯著升高。
圖12 磨刀門水道枯季表層和底層漲落潮平均鹽度分布Fig.12 Average salinity curves at water surface along Modaomen waterway
圖13 不同潮型下的縱剖面鹽度等值線Fig.13 Vertical section of salinity contours with different types of tides
3) 中潮漲急時, 河段1底層約7 km范圍內(nèi)由外海高鹽度水體占據(jù), 表層則由沖淡水覆蓋, 與洪灣水道匯合口位置的最大鹽度為5左右; 河段2在燈籠水閘以下游河段的鹽度變化范圍在 1~5, 底層最遠入侵距離達到燈籠水閘, 燈籠水閘往上游河段 2部分、河段3和河段4的水體鹽度都小于1。
4) 中潮落急時, 河段 1為沖淡水占據(jù), 呈現(xiàn)顯著的表、底分層現(xiàn)象, 底層最高鹽度達 17以上, 表層鹽度為 4左右, 與洪灣水道匯合口位置的鹽度達到6左右; 河段2的鹽度較漲急時整體有所升高, 底層 1的等值線向上游最遠接近聯(lián)石灣水閘, 從聯(lián)石灣水閘往上游河段2的部分、河段3和河段4的水體鹽度都小于1。
5) 小潮漲急時, 河段 1底層高鹽水體入侵較為顯著, 盡管底層含鹽度20以上的純海水僅僅進入到離口門 4 km的范圍, 但河段 1的底層基本都被 15以上的高濃度鹽水占據(jù), 表層為沖淡水占據(jù), 形成顯著分層現(xiàn)象, 與洪灣水道匯合口位置最高鹽度為8左右; 河段 2在聯(lián)石灣水閘下游河道的水體鹽度變化范圍在 1~10, 較高含鹽度出現(xiàn)在洪灣水道匯合口位置底層, 顯示是河段 1底層部分較高含鹽度的水體侵入所致; 聯(lián)石灣水閘以上河段的水體鹽度都在1以內(nèi);
6) 小潮落急時, 河段 1水體含鹽度較漲急時有顯著下降, 底層最大鹽度為 15左右, 表層最低為 3左右, 形成顯著表底分層現(xiàn)象, 底層較高濃度的含鹽水體由于地形阻礙無法完全退出口門, 鹽度等值線呈現(xiàn)由海向陸的傾斜分布, 與洪灣水道交匯處的最大鹽度為6左右; 河段2水體含鹽度較漲急時有所下降, 聯(lián)石灣水閘下游河道水體鹽度變化范圍為1~6, 從聯(lián)石灣水閘往上游河道含鹽度都在1以內(nèi)。
綜上可見, 近河口的磨刀門水道由于其特殊的水深地形和潮流動力特征, 決定了口門鹽度入侵的變化特性:
1) 枯季大潮時, 潮動力占絕對主導(dǎo)作用, 外海高鹽度的純海水從口門顯著入侵到磨刀門水道, 對磨刀門水道下泄徑流形成了顯著頂托作用, 致使大潮時的鹽度等值線總體呈現(xiàn)垂向分布; 中潮和小潮時, 潮動力減弱, 徑流作用凸顯, 密度較小的淡水或沖淡水由表層向下游流動, 底層高濃度的鹽水由于地形及補償流動原因向上游入侵, 形成顯著的表底分層現(xiàn)象。
2) 從口門至洪灣水道的河段 1, 由于水深地形呈現(xiàn)由海向陸傾斜下降的形態(tài), 導(dǎo)致底層容易為高濃度鹽水占據(jù), 且不易完全退出口門; 與洪灣水道的匯合口處地形的抬高對來自外海的底層高鹽水體入侵具有顯著的阻擋作用。
3) 大潮時的鹽水入侵距離較遠, 鹽度為 1的等值線能達到竹銀以上, 中潮和小潮時則基本都在聯(lián)石灣水閘以下。
4) 大潮和中潮期間, 落潮時的鹽水向上游入侵距離反而較漲潮時更遠, 根據(jù)對流態(tài)分析顯示, 潮汐動力強時, 漲潮動力強勁且枯季淡水徑流動力較弱, 致使上游下泄的淡水徑流聚集在磨刀門水道上游河段內(nèi)無法下泄, 導(dǎo)致上游水道含鹽度迅速降低;落潮時, 一方面聚集在水道內(nèi)的大量沖淡水由表層迅速下泄, 底層高鹽水體由于補償流動向上游入侵;另一方面磨刀門水道水深地形變化迅速, 造成表底層水流的強烈紊動, 進一步加劇了底層高含鹽水體往上游方向擴散。小潮期間, 由于整個水道內(nèi)水流流速很小, 流態(tài)平緩, 紊動較弱, 總體仍呈現(xiàn)漲潮時入侵距離大于落潮。因此, 枯季磨刀門水道鹽水入侵的主要影響因素為地形和潮動力。
為探討枯季磨刀門水道潮流和鹽水入侵特性,本文建立了磨刀門水道三維潮流和鹽度數(shù)學(xué)模型,對2009年枯季近半月的潮流和鹽水運動進行了模擬和分析, 結(jié)果顯示:
1) 枯季由于上游徑流量小, 磨刀門水道總體漲、落潮流速都不大, 表層總體漲潮平均流速都在0.5 m/s以內(nèi), 總體落潮平均流速在0.8 m/s以內(nèi); 底層總體漲落潮平均流速都在0.5 m/s以內(nèi);
2) 枯季磨刀門水道表、底層流速反向的時間段一般出現(xiàn)在經(jīng)過第 1個長時間的落潮流之后的第 2個較短落潮時間段內(nèi), 沿水深方向流向發(fā)生變化的位置隨潮動力減弱而有向水面移動的趨勢;
3) 從鹽度的平面分布來看, 磨刀門水道近口門河段(與洪灣水道匯合口以下)總體呈現(xiàn)漲潮時, 磨刀門水道東側(cè)含鹽度較西側(cè)大, 落潮時東側(cè)含鹽度較西側(cè)小的趨勢。竹排沙淺灘東側(cè)水道含鹽度總體較西側(cè)小, 大排沙淺灘東側(cè)石岐水道含鹽度則明顯較西側(cè)磨刀門水道大;
4) 潮汐動力較強時(大潮和中潮), 落潮時的鹽水向上游的入侵距離反而較漲潮時更遠, 潮汐動力弱(小潮)時, 整個水道內(nèi)水流流速很小, 流態(tài)平緩,紊動較弱, 總體仍呈現(xiàn)漲潮時入侵距離大于落潮時的情況。因此枯季磨刀門水道鹽水入侵特性取決于水道地形和潮動力。
[1]譚維炎.鹽水楔運動規(guī)律的研究述評[J].水科學(xué)進展,1994, 5(2): 149-159.
[2]徐峰俊, 朱士康, 王華.伶仃洋水動力環(huán)境分析及治理策略探討[J].人民珠江, 2004(1): 11-14, 27.
[3]劉杰斌, 包蕓, 黃宇銘.豐、枯水年磨刀門水道鹽水上溯運動規(guī)律對比[J].力學(xué)學(xué)報, 2010, 42(6):1098-1103.
[4]包蕓, 劉杰斌, 任杰, 等.磨刀門水道鹽水強烈上溯規(guī)律和動力機制研究[J].中國科學(xué): G輯: 物理學(xué)力學(xué)天文學(xué), 2009, 39(10): 1527-1534.
[5]聞平, 陳曉宏, 劉斌, 等.磨刀門水道咸潮入侵及其變異分析[J].水文, 2007, 27(3): 65-67.
[6]陳榮力, 劉誠, 高時友.磨刀門水道枯季咸潮上溯規(guī)律分析[J].水動力學(xué)研究與進展 A輯, 2011, 26(3):312-317.
[7]羅丹.磨刀門咸潮測驗的實踐[J].水利技術(shù)監(jiān)督, 2011,19(5): 21-23.
[8]韓志遠, 田向平, 劉峰.珠江磨刀門水道咸潮上溯加劇的原因[J].海洋學(xué)研究, 2010, 28(2): 52-59.
[9]賈良文, 吳超羽, 任杰, 等.珠江口磨刀門枯季水文特征及河口動力過程[J].水科學(xué)進展, 2006, 17(1):82-88.
[10]金臘華, 徐峰俊.河口及近海水質(zhì)模擬[M].北京: 化學(xué)工業(yè)出版社, 2007.
[11]陸仁強, 何璐珂.基于Delft3D模型的近海水環(huán)境質(zhì)量數(shù)值模擬研究[J].海洋環(huán)境科學(xué), 2012, 31(6):877-880.
[12]Lesseer G, Roelvink J A, Van Kester J A T M, et al.Development and validation of a three-dimensional morphological model[J].Coastal Engineering, 2004,51: 883-915.