EGU hourly FF10 processing: HAP calculation error

[callen52]

We have started processing EGU emissions with hourly FF10 inventories based on CEMs, in which pollutants other than NOX and SO2 are based on a new variable in the inventories called HOURACT. With this processing method, HAP emissions are not being calculated properly. I'm attaching a test case (scripts and inputs) to illustrate.

The test case includes two INVTABLEs, invtable_2017_NATA_CMAQ_15jun2023_v9.txt and invtable_2017_NATA_CMAQ_07may2024_nf_v13.txt. The v9 INVTABLE is our regular one, which includes some HAPs which have a factor of less than 1. For example, this HAP is in our ptegu inventory, and the factor 0.4406 is applied to all CAS=7738945 emissions from the annual inventory:
CHROMHEX 7738945 Y 0.4406 Y N tons/yr Chromium VI, PM

In the v13 INVTABLE, the factors for all HAPs - including 7738945 - are set to 1.

If the test case is run with the v9 INVTABLE, and then run a second time with the v13 INVTABLE, it is expected that CHROMHEX emissions will be different between the two runs, while all other pollutants should have identical emissions in both runs. However, what I see is that HAPs other than CHROMHEX end up with lower emissions when using the v9 INVTABLE compared to the v13 INVTABLE. My theory is that with this type of hourly processing, SMOKE is getting the factors from the INVTABLE mixed up and applying the wrong factors to the wrong HAPs.

SMOKE_ptegu_hourly_FF10_test_case_for_UNC.zip

[hnqtran]

Could not replicate the reported issue.

Descriptions on conducted test cases:

Case 1: Running Annual_point_onetime_ptegu_12US1_2021hb_cb6_21k_20240520091405.csh using invtable_2017_NATA_CMAQ_15jun2023_v9.txt (v9 INVTABLE)

Case 2: Running Annual_point_onetime_ptegu_12US1_2021hb_cb6_21k_20240520091405.csh using invtable_2017_NATA_CMAQ_07may2024_nf_v13.txt (v13 INVTABLE)

Case 3: Running Annual_point_onetime_ptegu_12US1_2021hb_cb6_21k_20240520091405.csh using a modified invtable_2017_NATA_CMAQ_15jun2023_v9.txt with doubling CHROMHEX CHROMHEX 7738945 Y 2.0 Y N tons/yr Chromium VI, PM (v9x2 INVTABLE)

Case 4: Running Annual_point_onetime_ptegu_12US1_2021hb_cb6_21k_20240520091405.csh using a modified invtable_2017_NATA_CMAQ_15jun2023_v9.txt with reducing CHROMHEX CHROMHEX 7738945 Y 0.0001 Y N tons/yr Chromium VI, PM (v9x0.0001 INVTABLE)

Case 5: Running Daily_point_summer_ptegu_12US1_2021hb_cb6_21k_20240520091645.csh using invtable_2017_NATA_CMAQ_15jun2023_v9.txt (v9 INVTABLE)

Case 6: Running Daily_point_summer_ptegu_12US1_2021hb_cb6_21k_20240520091645.csh using invtable_2017_NATA_CMAQ_07may2024_nf_v13.txt (v13 INVTABLE)

Comparison method:
Outputs (pnts*.ncf for annual script; pthour*.ncf for daily script)from smkinven in each case were compared using IOAPI m3diff tool for overall comparison of max, min and mean value. Additionally, the netcdf outputs were dumped into ascii files for line-by-line comparison

Findings:

Case 1 vs Case 2: CHROMHEX emission in Case 1 is smaller than in Case 2 (expected result); NOX and SO2 identical between two cases

Case 1 vs Case 3: CHROMHEX emission in Case 1 is smaller than in Case 3 (expected result); NOX and SO2 identical between two cases

Case 1 vs Case 4: CHROMHEX emission in Case 1 is higher than in Case 4 (expected result); NOX and SO2 identical between two cases

Case 5 vs. Case 6: CHROMHEX emission in Case 5 is smaller than in Case 6 (expected result); NOX and SO2 identical between two cases

I noticed that the run script provided by Christine call emf scripts smk_pt_annual_onetime_emf.csh and smk_pt_daily_emf.csh from smoke4.6. My tests were conducted with scripts from smoke5.0. It is unclear if the difference in version of emf script has any significant in the findings.

[callen52]

The smoke4.6 scripts are identical to the smoke5.0 scripts, so that should be fine.

When comparing the different cases, did you check pollutants other than CHROMHEX, NOX, and SO2? In my testing, NOX and SO2 are the same (as expected), but there are several HAPs which are different but should be the same, e.g. LEAD, ARSENIC, XYLENES, and ETHYLBENZ.

If you're seeing no differences, I wonder if the issue is somewhere downstream of Smkinven, perhaps in Temporal? If you were to go all the way through Smkmerge and generate model-ready emissions (which is what I had been doing in my testing), how do those compare?

[hnqtran]

I'm seeing LEAD, ARSENIC, XYLENES and ETHYLBENZ identical between v9 and v13. I have not looked pass smkinven since I though HAP processing happens in smkinven.
I will try subsequent SMOKE program and update with the results.

[hnqtran]

The issue was found to be caused by code line 892 - 908 of subroutine RDFF10PD in SMOKE/src/smkinven/rdff10pd.f which read CEM data after annual emission file was processed

              DO V = 1, NCEMPOL          ! loop for CEMS when HOURACT is processed

                IF( CEMPOL ) THEN
                    CDAT = EANAM( V )
                    COD  = INDEX1( CDAT, NIPPA, CEMPOLS )
                    IF( COD < 1 ) THEN
                        COD  = INDEX1( CDAT, NIPPA, EANAM )
                        CIDX = INDEX1( CDAT, NINVTBL, ITNAMA )
                    ELSE
                        CYCLE  ! skip CEMS polls
                    END IF
                END IF

C.............  Set conversion factor from Inventory Table. Default is
C               1., which is also what is used in all but a handful of
C               special toxics cases.
                CONVFAC = ITFACA( SCASIDX( UCASIDX( CIDX ) ) )

In the above code, if V is a CEMPOL (i.e., any emission inventory species in the annual emission file such as CHROMHEX, LEAD, XYLENES other than the NOX and SO2 in CEM file), the CDAT = EANAM(V) is the corresponding inventory pollutant name.

The problem is that in the INVTABLE, there are multiple CAS numbers/names that share a same CDAT inventory pollutant name. For example, in the invtable_2017_NATA_CMAQ_15jun2023_v9.txt, there are 25 CAS mapped to CHROMHEX (each with different conversion factor). Thus, as RDFF10PD uses CDAT to check for the corresponding pollutant index CIDX in the INVTABLE, it may get the wrong index. Consequently, as CIDX is then used to get conversion factor CONVFAC in the INVTABLE, wrong CIDX index leads to wrong CONVFAC value.

The above issue does not occur when smkinven processes annual emission input using RDEMSPD subroutine (src/smkinven/rdemspd.f) since CDAT is taken directly as pollutant CAS number/name (e.g., 7738945 instead of CHROMHEX) from raw emission input, based on which CIDX and CONVFAC are looked up correctly from the INVTABLE.

Another related issue is the double application of CONVFAC in line 965 - 971 of subroutine RDFF10PD

                        IF( CEMPOL ) THEN
                            EMISVA( HS,T ) = CONVFAC * EMIS(S,COD) * TDAT( D,H )  ! Store data in emissions
                            DYTOTA( HS,T ) = CONVFAC * EMIS(S,COD) * TOTAL
                        ELSE
                            EMISVA( HS,T ) = CONVFAC * TDAT( D,H )  ! Store data in emissions
                            DYTOTA( HS,T ) = CONVFAC * TOTAL
                        END IF

Here, if the pollutant is NOX or SO2 (CEMPOL = False), its annual emission is replaced by hourly emission from CEM file. For other pollutants (CEMPOL = True), its annual emission EMIS(S,COD) is multiplied with CONVFAC then hourly factor TDAT. The problem is that CONVFAC was already applied to EMIS(S,COD) in annual emission processing step.

[hnqtran]

Solution for the two issues above is to make CONVFAC = 1.0 in subroutine RDFF10PD for pollutants with CEMPOL = True

     IF (CEMPOL) THEN
          CONVFAC = 1.0
     ELSE
          CONVFAC = ITFACA( SCASIDX( UCASIDX( CIDX ) ) )
     END IF

[hnqtran]

Tagged B.H. to review this bug and its solution.

[bokhaeng]

@hnqtran I agree. It does look like the bug applies the split factor twice to EMIS during the hourly process. Instead of adding CONVFAC = 1.0 for CEMPOL, I think you could just remove CONVFAC from the following logic to make it simple.

                    IF( CEMPOL ) THEN
                        EMISVA( HS,T ) = EMIS(S,COD) * TDAT( D,H )  ! Store data in emissions
                        DYTOTA( HS,T ) = EMIS(S,COD) * TOTAL
                    ELSE
                        EMISVA( HS,T ) = CONVFAC * TDAT( D,H )  ! Store data in emissions
                        DYTOTA( HS,T ) = CONVFAC * TOTAL
                    END IF