AtChem / AtChem2

Atmospheric chemistry box-model for the MCM
MIT License
59 stars 23 forks source link

Non-integer Stoichiometry #513

Closed AlfredMayhew closed 9 months ago

AlfredMayhew commented 9 months ago

Hi, I've been working on running AtChem2 with the GEOS-Chem mechanism and have run into an issue that I've come across before, so I thought I'd mention it here in case it was functionality that we could look at adding to AtChem2.

Many lumped chemical mechanisms will include chemical reactions where a non-integer number of products are formed (or less-commonly, a non-integer number of reactants are consumed). This is not true for the MCM due to its semi-explicit nature. For example A + B --> 0.25C + 0.5D. AtChem2 does not currently support writing reaction stoichiometry in this way, instead relying on duplicating species within the equation (e.g. A + B --> 2C is written as A + B --> C + C). Again, this works fine for the MCM but I believe it would add a lot more flexibility to AtChem2 to be able to accept numbers before species names as stoichiometry. A search of previous issues brought up this discussion from 2018 where you agreed that duplicating species was sufficient since it would work for the MCM.

I have got around this issue a couple of different ways in the past by breaking the reactions apart and duplicating the number of species as required. I think the most generalizable way I have found is to split the reaction based on each product, and then add an extra dummy reaction and products to account for the loss of the reactants. For example, % k : A + B --> 0.25C + 0.5D could be written as:

% k*0.25 : A + B --> C + A + B ;
% k*0.5 : A + B --> D + A + B ;
% k : A + B --> ;

This isn't a bad work-around but it does increase the size of the mechanism considerably if there are multiple lumped reactions, so it would be nice to be able to have AtChem2 read the stoichiometry of reactions. This approach would also mess up the rate files output from AtChem2. I'm happy to help look into how to implement this if you think it would be doable, I just have no idea where to start looking!

Thanks, Alfie

rs028 commented 9 months ago

Hi Alfie, this is a good point and I agree it would be useful to have that flexibility. It would also perhaps avoid a possible source of error if the .fac file is not formatted properly.

One way to do it is to pre-process the fac file and split the chemical equation, as you suggest. Or it can be done by adding stoichiometric coefficients to the Fortran subroutine. The second option is more robust perhaps.

I think we can be fairly certain that no chemical species will have a name beginning with a number and therefore if there is an initial number it can be safely assumed that it corresponds to a stoichiometric coefficient?

AlfredMayhew commented 9 months ago

Yes, I agree that it's reasonable to say that species names will not begin with a number. This has been true for all mechanism formats and modelling frameworks I've come across.

I think pre-processing the fac file would mean we'd have to adjust the rate outputting so that the split reactions aren't output. It seems less convoluted to carry the stoichiometry directly through to the calculations, but I guess it could be a bigger job. I'll have a dig around if I get some time and let you know if I make any progress.