## Concentration / Display / Source-Receptor Solve

Overview: This menu can be used to solve the source-receptor coefficient matrix for the source term vector given a measured data vector and where the matrix values are the dilution factors for each source-receptor pair as created through the Run Matrix menu option. This requires that the ICHEM=1 option is set in the SETUP.CFG namelist file. The measured data vector at multiple receptor locations is required and it must be defined in the DATEM format. The format of this file is discussed more detail in the GeoLocation menu. Due to model errors and insufficient sampling not all simulations will yield a solution. Results are written to the file: matrix.txt

Theory: Assume that the concentration at receptor R is the linear sum of all the contributing sources S times the dilution factor D between S and R:

- S1D11 + S2D12 = R1
- S1D21 + S2D22 = R2

The dilution factors are defined as the coefficient matrix. A plot of the column product SiDij can be presented as a map of the concentrations contributed by source i to all the receptors. A plot of the row product SiDij for receptor j would show a map of the concentration contributed by each source to that receptor. In this situation it is assumed that S is known for all sources. The elements of the coefficient matrix are normally computed explicitly from each source to all receptors, the traditional forward downwind calculation.

In the case where measurements are available at receptor R and source S is the unknown quantity, the linear relationship between sources and receptors can be expressed by the relationship:

- Dij Si = Rj,

which has a solution if we can compute the inverse of the coefficient matrix:

- Si = (Dij)-1 Rj.

For example, in the case of an unrealistic 2x2 square matrix (the number of sources equals the number of receptors), the inverse of D is given by:
| +D22 -D12 | 1/(D11D22-D12D21)
| -D21 +D11 |

The solution for S1 (first row) can be written:

- D22R1/(D11D22-D12D21) - D12R2/(D11D22-D12D21)

As a further simplification, assume that there is no transport between S1 and R2(D12 = 0), and then the result is the trivial solution that the emission value is just the measured concentration divided by the dilution factor:

- S1 = R1/D11

The matrix solution has three possibilities. The most common one is that there are too many sources and too few receptors which results in multiple solutions requiring singular value decomposition methods to obtain a solution. The opposite problem is that there are too many receptors and too few unknowns hence an over determined system requiring regression methods to reduce the number of equations. Unfortunately, possibilities for a matrix solution may be limited at times due to various singularities, such as columns with no contribution to any receptor or measured values that have no contribution from any source. The solution to these problems is not always entirely numerical as the underlying modeling or even the measurements can contain errors. Note that large dilution factors (very small predicted concentrations) at locations with large measured values will lead to large emissions to enable the model prediction to match those measurements. The opposite problem also exists in that negative emission values may be required to balance high predictions with small measurements. The solution to the coefficient matrix is driven by errors, either in the measurements, the dispersion model, or the meteorological data.

Step 1: defines the binary matrix file which is the output file from the matrix simulation that produced model predicted concentrations at an output interval that can be matched to the measured sampling data following the procedure described in the Run Matrix menu option. The resulting coefficient matrix should be viewable through the Display/Matrix menu tab. Ideally the model simulation emission rate should be set to a unit value. The input examples\matrix\control_srm file and examples\matrix\setup_srm file should be retrieved into the setup and configuration menus to create the binary file. This CONTROL file has been pre-configured for a regular matrix of 81 sources between 38N-42N and 82W-78W at 0.5 degree resolution. The model can be run from the Run Standard Model menu tab when using this file. If the 3-source location method is used to configure the simulation, then it should be run from the Special Runs / Matrix menu tab. Also note that there is no requirement that the source vector (or the sampling vector) of the coefficient matrix be defined on a regular latitude-longitude grid.

Step 2: defines the measured data input file which is an ASCII text file in the DATEM format. The first record provides some general information, the second record identifies all the variables that then follow on all subsequent records. There is one data record for each sample collected. All times are given in UTC. This file defines the receptor data vector for the matrix solution. It may be necessary to edit this file to remove sampling data that are not required or edit the simulation that produces the coefficient matrix to insure that each receptor is within the modeling domain.

Step 3: defines the receptor to source unit conversion from sampling units to the emission units. For instance, the default value of 10-15 converts sampling data measured in pg/m3 (pico = 10-12) to emission units of kg. There are 10+15 pg in one kg.

Step 4: solves for the source vector using the SVD (Singular Value Decomposition) methods from Numerical Recipes, The Art of Scientific Computing, Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1986, Cambridge University Press, 818 p. There are no options that control the solution, but the verbose checkbox shows the coefficient matrix during each iteration of the solution in addition to the final source vector results. Each iteration results in additional automatic editing of the coefficient matrix to remove or reduce singularities. The iteration process will continue until the number of unknowns (sources) equals the number of equations (samplers) or until the iteration limit is reached. With each iteration, singular columns are removed, and the lowest coefficient values are also eliminated at factor-of-two intervals from the lowest value. In the event the computation fails to converge on a solution, it may be necessary to either edit the measured data file to limit the observations or modify the control file to change the location or resolution of the source grid.

Step 5: displays the solution results from the output file: matrix.txt. The output detail depends upon the verbose checkbox and the number of iterations required. A solution may contain negative values as well as extreme positive emission results. Such values are not realistic and are a result of model errors or other uncertainties discussed previously. Only positive values within 6 orders-of-magnitude from the minimum emission value are shown. Below are the results from the 2nd pass in the output file matrix.txt for the example calculation. Note that due to precision issues, processor type, and operating system, the actual numbers you may get will vary. The values shown below were solved on a 64 bit operating system.

- EQN = VARB Source vector solution pass #2
- 41.0 -82.0 0.0
- 40.5 -81.5 0.0
- 41.0 -81.5 0.0
- 40.5 -81.0 4.482
- 40.0 -80.5 0.0
- 40.5 -80.5 87.321
- 40.0 -80.0 6.152
- 40.5 -80.0 0.0
- 39.5 -79.5 0.0
- 40.0 -79.5 0.0
- 39.5 -79.0 2.515
- 40.0 -79.0 0.0
- 39.5 -78.5 0.0
- 39.0 -78.0 0.0
- 39.5 -78.0 0.0
- Average: 6.698

It is not surprising that multiple locations show realistic positive emissions due in part to the fact that multiple source locations contribute (non-zero coefficient matrix values) to many of the 6-hour duration samples. However as can be seen in the illustration below, when the model estimated emissions are represented spatially, their alignment is along the primary transport axis (NW to SE) with a spatial extent comparable to the high probability region illustrated in the GeoLocation calculation for this same example case. Clearly the modeling and sampling protocol is insufficient to determine the exact location and emission amount (10 kg) with any more precision, but the emission value averaged over the 15 locations is 6.7 kg.

Alternate computational approaches might be considered by computing only those elements of the matrix for which receptors that have a non-zero value (and perhaps a few surrounding zero points), using the upwind dispersion calculation approach. Although the advection term of the dispersion equation is reversible, the diffusion term is not, and therefore these calculations produce a field that looks like concentration and of similar magnitudes. In this respect, the source attribution function could be treated in a similar manner as the dilution factor elements in the coefficient matrix. In this case the dilution factor is not the dilution of the source to the receptor, but the dilution of the upwind source regions to the receptor under consideration. One positive outcome of such an approach is that each receptor point is assured of contributing a dilution factor to one or more source locations and eliminating the possibility of singularities. Zero measurements at sampling locations outside of the measured plume or those that occur before or after the plume passage can be added to the matrix as long as they contribute positive dilution factors to the source vector. Although not yet available through the GUI, such an option may be added to a future version.

Table of Contents