增量濃度分析程序與檢討

前言/背景

計算不同排放情境模擬結果的濃度差異，是執行空品模式常見的作業。如CMAQ這樣有nc檔案輸出的模式，以dNC就可以簡單解決。但因為CMAQ的粒狀物定義同時有i,j,k濃度(CCTM_ACONC)與其重量比例(CCTM_APMDIAG)等2個場同時介入，因此造成非常高度之非線性結果。
按照PM的計算方式(see also [[2021-12-16-combine]][^1])，非線性來源有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在某一時間之山區、電廠下游的西南海域遠方出現NO₃^-增量、而在海面則出現負值之增量(空白處為無法取log值之負值增量區域)，為化學反應造成的非線性現象，即使取月均值，電廠附近局部網格仍有負值。其他SO₄⁼也出現負值增量的結果。
營運前後粗細粒比例(CCTM_APMDIAG檔)具有差異性，如圖2a以營運前為1.0，營運後PM25AC月均值的增加幅度，圖2b則為PM10AC的增加幅度。因電廠運轉造成煙流範圍細粒比例增加約1%，洽與上圖₃^-負值增量及PM_2.5與PM₁₀增量之差值負值分布相同
PM_2.5增量高於PM₁₀的增量(圖3)，圖中只顯示負值部分。顯示電廠為起點之煙流範圍有負值之最大值，約與上圖NO₃^-負值增量濃度分布相同，為combine的最終結果，包括污染項目間的集成與PM分率的效果。


圖1a 2019/01/01/00Z 興達新3氣機組PM_2.5中NO₃^-濃度值之增量。	圖1b 同左但為月平均值

圖2a 計畫營運前後Aitken mode濃度在PM_2.5部分之比例(PM25AC)之月均值增加率。	圖2b 同左但為PM₁₀部分比例(PM10AC)。電廠造成附近PM10AC月均值增加0.01%

圖3a 同圖1a時間PM_2.5與PM₁₀增量之差值	圖3b 同左但為月平均值

檢討

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

出現負值增量與PM₁₀小於PM_2.5可能的原因列表

項目	公版模式指定內容及計算公式	出現負值增量與PM₁₀小於PM_2.5可能的原因	此次修正作法
合併程序	combine.sh	公版只提供combine後之背景濃度，未提供base之ACONC及APMDIAG檔，因此只能先進行combine，結果再進行差異分析	先進行營運前後、個別成分逐時濃度之篩選、再進行combine與差異分析
物種定義檔，執行combine讀取	SpecDef_cb6r3_ae7_aq.txt		(未修正)
PM₁₀物種組成定義於物種定義檔	ATOTI[0]×PM10AT[3] + ATOTJ[0]×PM10AC[3] + ATOTK[0]×PM10CO[3]	負值可能來自相乘的元素及加總過程	(未修正)
PM25_TOT物種組成定義於物種定義檔	ATOTI[0]×PM25AT[3] + ATOTJ[0]×PM25AC[3] + ATOTK[0]×PM25CO[3]	負值可能來自相乘的元素及加總過程	(未修正)
重做base之小時平均濃度輸出檔CCTM_ACONC	提供前述公式中各I,I,K各mode成份之總合ATOTI, ATOTJ, ATOTK	因化學反應case營運後部分逐時濃度結果反而低於base個案，PM₁₀成分累積更多負值造成低於PM_2.5之結果	過濾營運後逐時濃度低於背景之結果(設為背景值)。
網格點上各粒徑比例之小時平均值CCTM_APMDIAG	提供前述公式中之PM10AT, PM10AC, PM10CO，以及PM25AT, PM25AC, PM25CO	營運前後部分地區時間的粒徑比例發生差異，增量煙流之細粒比例有所增加。如該處出現負值增量，則會對該處平均值放大扣減效果。	只選擇營運後之粒徑比例進行combine計算，避免不一致。

因應策略方案

先相減再過濾

雖然就化學反應而論，負值增量是合理與必然的結果，但是在法規應用上為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/* .

計算程序比較

順序	傳統作法	新計算程序建議	比較說明
1	營運前後模擬結果分別進行combine	營運前後CCTM_ACONC相減、過濾^*	後者檔案較大、需要較大磁碟空間進行作業、妥善檔案管理、平行運作
2	營運前後CCTM_COMBINE檔案再進一步篩減項目(shk.cs)	連結營運後之CCTM_APMDIAG檔案到同一目錄	前者可以再進一步縮小檔案容量
3	營運前後shk.cs結果相減、分析、繪圖	營運前後以同一組CCTM_APMDIAG進行combine	前者完工
4	(無)	進行shk.cs、分析、繪圖	後者完工

^*：先相減後過濾(此處作法)、或為套用公版後處理作業架構先過濾再進行combine、相減，效果應為一致。
結果及討論
PM_2.5：最大值變動不大，但臺灣本島中北部、東部的增量變多了。海面上、電廠附近的負值增量消失了。
PM₁₀：山區低濃度範圍更擴大一些
PM_2.5/PM₁₀：月均值的比例，未過濾前的振盪很大、包括在外海、北部地區及山區，電廠附近也有負值之振盪。經過濾後模擬範圍的比例在0.35~0.95之間，山區約在0.7上下，平地範圍較高約在0.8~0.9，越遠越高。電廠附近與煙流也較高。