增量濃度分析程序與檢討

前言/背景

• 計算不同排放情境模擬結果的濃度差異，是執行空品模式常見的作業。如CMAQ這樣有nc檔案輸出的模式，以dNC就可以簡單解決。但因為CMAQ的粒狀物定義同時有i,j,k濃度(CCTM_ACONC)與其重量比例(CCTM_APMDIAG)等2個場同時介入，因此造成非常高度之非線性結果。
• 按照PM的計算方式，非線性來源有3：
  1. 化學反應：開發單位(以下以高雄興達電廠為例)因排放大量NO，造成局部地區O₃快速反應消失，因氧化物減少，原生性污染物SO₂、NO₂的轉化能力下降，於是造成局部地區衍生性粒狀物成分減量、造成PM比例的失衡，會出現在營運前後CCTM_ACONC檔案的差異。
  2. 營運前後PM粒徑重量比例係數(CCTM_APMDIAG)的差異。營運後不必然(也不可能)所有時間、所有地區、所有mode的比例都是增加。
  3. PM集成的過程：不論是PM_2.5或PM₁₀的計算過程中，都會先分別針對i,j,k mode將所有粒狀物進行加總。如果有的物質增加、有的減少，總合後將會使營運前後的差值不如預期(不必然增加)。
• 由於目前還沒有找到討論CMAQ電廠個案PM_2.5/PM₁₀增量比例的期刊文章，整體計算方式值得詳細討論。

傳統作法

• 由於CCTM_ACONC檔案的化學物質多達200餘種，再加上CCTM_APMDIAG逐日結果，檔案非常龐大，需先進行整併以減省空間。
• 因此大多數處理流程會先進行個案模擬的combine，再進行相減(dNC)與增量的後續統計。

結果

1. 相減結果會在看似非關連之區域出現微幅的濃度擾動。如圖1在某一時間之山區、電廠下游的西南海域遠方出現NO₃^-增量、而在海面則出現負值之增量(空白處為無法取log值之負值增量區域)，為化學反應造成的非線性現象，即使取月均值，電廠附近局部網格仍有負值。其他SO₄⁼也出現負值增量的結果。
2. 營運前後粗細粒比例(CCTM_APMDIAG檔)具有差異性，如圖2a以營運前為1.0，營運後PM25AC月均值的增加幅度，圖2b則為PM10AC的增加幅度。因電廠運轉造成煙流範圍細粒比例增加約1%，洽與上圖₃^-負值增量及PM_2.5與PM₁₀增量之差值負值分布相同
3. PM_2.5增量高於PM₁₀的增量(圖3)，圖中只顯示負值部分。顯示電廠為起點之煙流範圍有負值之最大值，約與上圖₃^-負值增量濃度分布相同，為combine的最終結果，包括污染項目間的集成與PM分率的效果。

圖1a 2019/01/01/00Z 興達新3氣機組PM2.5中NO3-濃度值之增量。	圖1b 同左但為月平均值
圖2a 計畫營運前後Aitken mode濃度在PM2.5部分之比例(PM25AC)之月均值增加率。	圖2b 同左但為PM10部分比例(PM10AC)。電廠造成附近PM10AC月均值增加0.01%
圖3a 同圖1a時間PM2.5與PM10增量之差值	圖3b 同左但為月平均值

檢討

• 影響PM結果的因素檢討起來有下列3項：
  • 增量中含有之負值、且/或
  • PM25AC、PM10AC之增加、且/或
  • combine的集成效果
• 3因素的相依性檢討
  • 即使使用同一組CCTM_APDIAG（不論營運前或營運後），因一般NO₃^-在粗顆粒佔比較大，負值則將造成粗顆粒的減少，因此造成PM₁₀增量小於PM_2.5。
  • 即便ATOT_i為正值，營運後乘上較大的PM25AC比例得到較高的PM_2.5，是有可能大於乘上增加比例不太多的PM10AC的結果。
  • 由於ATOT的計算公式是相加，唯有每個ATOT成分都滿足營運後大或等於營運前，才能保證一定能得到ATOT也是營運後大或等於營運前的結果。
• 此3者同時作用
  • PM_2.5增量大於PM₁₀雖不合理，但其乃為必然之數學結果。

因應策略方案

先相減再過濾

• 雖然就化學反應而論，負值增量是合理與必然的結果，但是在法規應用上為trivial solution，無法討論。
• 且在時間平均過程中負值會造成干擾、降低平均結果而不保守。
• 建議在combine之前就將負值排除，以避免粒徑分率與集成過程的干擾。

#kuang@DEVP /nas2/cmaqruns/2019force/output/2019-01/grid03/cctm.XindaN3G-BASE_withFilter0
d1=../../cctm.BASE/daily
d2=../../cctm.XindaN3G/daily
for i in {01..31};do dNC $d1/CCTM_ACONC$i.nc $d2/CCTM_ACONC$i.nc CCTM_ACONC$i.nc ;done

• 使用簡單的np.where即可完成過濾的動作。

V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
for v in V[3]:
  var=np.where(nc[v][:]>0,nc[v][:],0)
  nc[v][:]=var[:]

先過濾再相減

• 如果為符合公版模式現行作業程序(需同時提供營運前後combine結果)，要在相減前要進行過濾
• 設若nc為營運後、nc0為營運前之背景濃度，前述np.where條件需改成「營運後濃度是否大於背景」，如否，則設成背景濃度。

V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
for v in V[3]:
  var=np.where(nc[v][:]>nc0[v][:],nc[v][:],nc0[v][:])
  nc[v][:]=var[:]

必須使用同一組粒徑分率(CCTM_APMDIAG)進行combine

• 如果在combine前先進行濃度差值的計算，增加的濃度應該適用營運後的CCTM_APMDIAG比較合理。

#kuang@DEVP /nas2/cmaqruns/2019force/output/2019-01/grid03/cctm.XindaN3G-BASE_withFilter0
ln -s ../../cctm.XindaN3G/daily/APDIAG_b_n3g/* .

計算程序比較

• ^*：先相減後過濾(此處作法)、或為套用公版後處理作業架構先過濾再進行combine、相減，效果應為一致。

  結果及討論

• PM_2.5：最大值變動不大，但臺灣本島中北部、東部的增量變多了。海面上、電廠附近的負值增量消失了。
• PM₁₀：山區低濃度範圍更擴大一些
• PM_2.5/PM₁₀：月均值的比例，未過濾前的振盪很大、包括在外海、北部地區及山區，電廠附近也有負值之振盪。經過濾後模擬範圍的比例在0.35~0.95之間，山區約在0.7上下，平地範圍較高約在0.8~0.9，越遠越高。電廠附近與煙流也較高。

圖4a 2019年1月份平均PM2.5增量濃度	圖4b 同左但為新計算程序
圖5a 2019年1月份平均PM10增量濃度	圖5b 同左但為新計算程序
圖6a 2019年1月份平均PM2.5增量濃度與PM10之比值	圖6b 同左但為新計算程序