This program is created to integrate X-STILT and OCO-3 SAM observations in the Bayesian inversion to optimize the point source CO2 emissions (ex. the power plant).
git clone --recursive https://github.com/dogtorXiao/WRF-XSTILT_Inversion.gitIn order to run the demo_run (the main data and figures shown in the paper "High-Resolution Modeling to Quantify CO2 Emissions from Industrial Point Sources"), please set a virtual environment with R packages by:
install.packages(c("reshape2", "ggplot2", "dplyr", "raster", "ggmap", "lubridate", "patchwork"))
in R session.
Enter the file demo_run/xstilt_post_dealing.r, set the case index (int) and the spriorisf (float) (the ratio between the assumed prior uncertainty and the prior emission), insert your google api key.
To run, in a terminal:
cd demo_run
Rscript xstilt_post_dealing.r
cd InvResults
Rscript plot.chart_converted.r
The output figures can be found in demo_run/DataFigures/YYYYMMDD
- Recommended: conda environment with python3.11 and package requirements in requirements.txt
- Interactive GUI session is in need for plume segmentation
- R interpreter, prerequisited packages will be automatically installed by submodule X-STILT, make sure the python interpreter can be called by R system() function
- meteorological data in arl format (eg. https://www.ready.noaa.gov/data/archives/gfs0p25/, or convert the netcdf meteorological data files to arl format (https://www.ready.noaa.gov/HYSPLIT_data2arl.php#INFO))
- radiosonde data (https://ruc.noaa.gov/raobs/)
- ODIAC emission (https://db.cger.nies.go.jp/dataset/ODIAC/DL_odiac2020b.html)
- OCO-3 CO2 observations (https://search.earthdata.nasa.gov/search?q=OCO3_L2_Lite_FP)
- TROPOMI NO2 data (https://search.earthdata.nasa.gov/search?q=S5P_L2__NO2____1)
STEPS:
Shared parameters setting:
key: Input a Google Map API key and insert into filesite,site_lon_lat: Find the geolocation information in Google Maps (https://www.google.com/maps) - find the latitude and longitude, then setsiteandsite_lon_lattimestr: OCO-3 overpass time in 'YYYYMMDD'trp_path: TROPOMI data pathraob_path: Radiosonde data pathtiff.path: ODIAC data pathmet_path: ARL format meteorology data path
-
Use the X-STILT forward function to get a sense of how the downwind domain looks like (usually just take a few minutes)
run_forward = T:box.len,dtime.from,dtime.to,dtime.sep,nhrs_forwardrequired
-
Plot the downwind domain, segment the inside pixels, and click the pixels farthest away from the point sources (up to 40) as a set of samples. Get the X-STILT backgrounds and background errors.
run_bg_xstilt = T:tdrequired
-
Segment the forward plume downwind domain (manually click the plume edge, interactive GUI needed):
run_seg = T
-
Release the particles at the samples to run for plenty of time with X-STILT backward function
run_sample_xstilt = T
-
Calculate the time period needed to run X-STILT and the average PBL height. Determine the scaling factor for the PBL height (zisf) with any observations you have (if there is no data available, skip this step)
run_PBL_nhrs_calc = T
MAIN_RUN = T
-
Decide the standard zisf and turbulented zisfs to run X-STILT for multiple times (vertical error)
run_ver_main = T:zisfs(vector, the standard and turbulent scaling factors of the PBL heights)zisf_indx(the index of zisfs)nhrs_backminaglmaxagl
-
The wind error running (using the radiosonde data) (transport error)
run_wind_main = T
POST_RUN = T
- Calculate the simulated XCO2 and transport errors (With ODIAC emission inventory as the prior), merge the observed, convolved, and turbulented XCO2, plot the XCO2 and observation minus simulation difference
bg: the CO2 background
INVERSION_RUN = T
Parameters setting:
bg,bg_uncert: the CO2 background and uncertainty obtained from NO2 mask or X-STILT forward functionspriorisf: the ratio between the prior uncertainty and the prior emissionslength_scale_priori: prior uncertainty length scale in kmlength_scale_obs: observation uncertainty length scale in kmdomain_opt_deg: box area for the optimization as the prior emission vector (in degree)
-
Extract the NO2 data within the
bg_radius(in km), select the NO2 plume, manually click the NO2 plume (interactive GUI required), extract the OCO-3 CO2 pixels inside the NO2 plumerun_no2_dealing = T:buff_deg,bg_radiusfor NO2 plume masking
-
Conduct the inversion (The point source is very sensitive compared to area fluxes, while the SAM mode observations provides more than 1000 pixels around the point source, most of which are easy to introduce uncertainties. Therefore, selecting the pixels in the inversion is crucial. Here we provide 2 methods:
Method 1: Only select the pixels inside the near-time TROPOMI NO2 plumes. Before conducting the inversion, please plot the NO2 plumes and segment the OCO-3 pixels, calculate the background errors
Method 2: Select all the pixels in the downwind domain)
run_inversion = T:inversion_recp = recp_fn: inversion using all pixels in the X-STILT forward downwind domaininversion_recp = recp_fn_in_NO2: inversion using only pixels in the NO2 plume mask