Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

H2 effects on CH4 lifetime #758

Open
wants to merge 7 commits into
base: dev-h2
Choose a base branch
from
+138 −74
Open

H2 effects on CH4 lifetime #758

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2203297
use more conservative code to define a constant
kdorheim Nov 4, 2024
a46e14e
increase the version index so that it is easier to keep track of erro…
kdorheim Nov 6, 2024
ec1511f
implement the H2 coeff for the OH componet as something that can be r…
kdorheim Nov 6, 2024
42479b7
saving wip
kdorheim Nov 11, 2024
8ed3f3e
Correct warning message
kdorheim Nov 26, 2024
95b3196
feedback from ben
kdorheim Dec 2, 2024
7058e56
address response to pr 758
kdorheim Dec 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Open in desktop
kdorheim marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: hector
Title: The Hector Simple Climate Model
Version: 3.2.0
Version: 3.3.0.9999999
kdorheim marked this conversation as resolved.
Show resolved Hide resolved
Authors@R: c(person("Kalyn", "Dorheim",
email = "[email protected]",
role = c("aut", "cre"),
Expand Down
4 changes: 4 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
# hector 3.3.0.9999999

* under development for the H2 integration, increasing the version index though will be ensure that the right version of hector is being run with gcam

# hector 3.2.0
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10698028.svg)](https://doi.org/10.5281/zenodo.10698028)
* Correct aerosol forcing coefficients based on Zelinka et al. (2023)
Expand Down
7 changes: 7 additions & 0 deletions R/RcppExports.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@
#' @describeIn emissions Emissions hydrogen gas (\code{"Tg H2"})
NULL

#' @describeIn oh coefficent for the h2 emissions
NULL

#' @describeIn msgtype Message type for retrieving data from a component
#' @keywords internal
GETDATA <- function() {
Expand Down Expand Up @@ -1571,6 +1574,10 @@ HEAT_FLUX <- function() {
.Call('_hector_HEAT_FLUX', PACKAGE = 'hector')
}

COEFF_H2 <- function() {
.Call('_hector_COEFF_H2', PACKAGE = 'hector')
}

newcore_impl <- function(inifile, loglevel, suppresslogging, name) {
.Call('_hector_newcore_impl', PACKAGE = 'hector', inifile, loglevel, suppresslogging, name)
}
Expand Down
1 change: 1 addition & 0 deletions inst/include/component_data.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -281,6 +281,7 @@
#define D_COEFFICENT_CH4 "CCH4"
#define D_COEFFICENT_NMVOC "CNMVOC"
#define D_COEFFICENT_CO "CCO"
#define D_COEFFICENT_H2 "CH2"
#define D_EMISSIONS_H2 "H2_emissions"


Expand Down
2 changes: 1 addition & 1 deletion inst/include/h_util.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
* \brief The model version number to be included in logs and outputs.
* \note Manually update the git tag to match this.
*/
#define MODEL_VERSION "3.2.0"
#define MODEL_VERSION "3.3.0.9999999"

#define OUTPUT_DIRECTORY "output/"

Expand Down
2 changes: 2 additions & 0 deletions inst/include/oh_component.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ class OHComponent : public IModelComponent {
double CNMVOC; // coefficent for NMVOC
double CNOX; // coefficent for NOX
double CCH4; // coefficent for CH4
unitval CH2; // coefficent for CH4
kdorheim marked this conversation as resolved.
Show resolved Hide resolved


// logger
Logger logger;
Expand Down
1 change: 1 addition & 0 deletions inst/input/hector_ssp119.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same spelling note; also, is there a reference/source for this value? If so please note it.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was there some derivation needed for this, or was this coefficient in the right format from the original source? (And how similar were the coefficients for other emissions to the ones we are using now?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay there was a derivation that was needed for this parameter, which will be described in detail in the manuscript associated with this work although that is a TBD, so I have included a TODO-H2 tag which is linked to the H2 project as a reminder to address this issue before merging into main or anything.

@ssmithClimate Could you elaborate on this a bit more? (And how similar were the coefficients for other emissions to the ones we are using now?) I am not sure if I understand what you are asking here. It is also unclear where the original coefficients came from @ssmithClimate do you know the source of those?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the delay. The coefficients in hector now ultimately derive from this paper:
Wigley, T.M.L., Steven J. Smith, and M.J. Prather (2002) Radiative Forcing due to Reactive Gas Emissions. Journal of Climate 15(18), pp. 2690–2696

Although the hector paper quotes the source as :
Tanaka, K., Kriegler, E., Bruckner, T., Hooss, C., Knorr, W., and Raddatz, T.: Aggregated Carbon Cycle, Atmospheric Chemistry, and Climate Model (ACC2) – description of the forward and inverse models, 1–188, Max Planck Institute for Meteorology, Hamburg, Germany, 188, 2007.
(which is a more indirect reference)

So my question was - is the overall chemistry representation of the source we're drawing the H2 coefficients from is generally consistent with the chemistry representation used now in Hector. If not it's something to flag (e.g., means we ultimately might need to update more than H2) since the two things we're merging together are not entirely consistent with each other. (Given uncertainty its likely fine, but we would want to flag this at least.)

So the way I thought to evaluate that was to look if the relative magnitude of any common coefficients for the non-H2 components of the equations that were used in the H2 work we're using now and the original equation we use now in Hector were similar. For example is the ratio of the impact of NOx to impact of CH4 or NOx to CO in the original hector (e.g. Wigley et al, Tanaka et al) are similar to what's in the H2 papers we are citing (to the extent they contain these other effects).

For example, in hector we have something like this for the CH4 OH lifetime (Hartin et al. 2015):

ln(OH)t = −0.32(ln[CH4 ]t − ln[CH4 ]t 0 ) + 0.0042(E(NOx)t ) − (E(NOx)t0) − 0.000105(E(CO)t ) − (E(CO)t0) − 0.00315(E(VOC)t ) − (E(VOC)t0)

So the impact of NOx vs CO is: 0.0042 / 0.000105 = 40 times (on an emissions basis in this case). Or we could look at the relative impact of NOx to CH4 in this equation if that's what's common between what's in hector now and what is represented in the equations from the paper's we are drawing from for the new H2 components.

To do this we need to have both sets of equations cast into an equivalent form.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kdorheim comments here!


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp126.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp245.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp370.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp434.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp460.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp534-over.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
2 changes: 2 additions & 0 deletions inst/input/hector_ssp585.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,8 @@ CNOX=0.0042 ; coefficent for NOX
CCO=-0.000105 ; coefficent for CO
CNMVOC=-0.000315 ; coefficent for NMVOC (non methane VOC)
CCH4=-0.32 ; coefficent for CH4
CH2=-0.00044625 ; coefficent for H2


;------------------------------------------------------------------------
[ozone]
Expand Down
11 changes: 11 additions & 0 deletions src/RcppExports.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2627,6 +2627,16 @@ BEGIN_RCPP
return rcpp_result_gen;
END_RCPP
}
// COEFF_H2
String COEFF_H2();
RcppExport SEXP _hector_COEFF_H2() {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
rcpp_result_gen = Rcpp::wrap(COEFF_H2());
return rcpp_result_gen;
END_RCPP
}
// newcore_impl
Environment newcore_impl(String inifile, int loglevel, bool suppresslogging, String name);
RcppExport SEXP _hector_newcore_impl(SEXP inifileSEXP, SEXP loglevelSEXP, SEXP suppressloggingSEXP, SEXP nameSEXP) {
Expand Down Expand Up @@ -3035,6 +3045,7 @@ static const R_CallMethodDef CallEntries[] = {
{"_hector_FLUX_MIXED", (DL_FUNC) &_hector_FLUX_MIXED, 0},
{"_hector_FLUX_INTERIOR", (DL_FUNC) &_hector_FLUX_INTERIOR, 0},
{"_hector_HEAT_FLUX", (DL_FUNC) &_hector_HEAT_FLUX, 0},
{"_hector_COEFF_H2", (DL_FUNC) &_hector_COEFF_H2, 0},
{"_hector_newcore_impl", (DL_FUNC) &_hector_newcore_impl, 4},
{"_hector_shutdown", (DL_FUNC) &_hector_shutdown, 1},
{"_hector_reset", (DL_FUNC) &_hector_reset, 2},
Expand Down
2 changes: 1 addition & 1 deletion src/ch4_component.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ void CH4Component::run(const double runToDate) {
<< std::endl;

// Permafrost thaw produces CH4 emissions
#define PG_C_TO_TG_CH4 (1000.0 * 16.04 / 12.01)
#define PG_C_TO_TG_CH4 (1000.0 * 16.04 / 12.01);
const double rh_ch4 =
core->sendMessage(M_GETDATA, D_RH_CH4).value(U_PGC_YR) * PG_C_TO_TG_CH4;

Expand Down
27 changes: 22 additions & 5 deletions src/oh_component.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,12 +57,15 @@ void OHComponent::init(Core *coreptr) {

// Inform core what data we can provide
core->registerCapability(D_LIFETIME_OH, getComponentName());
core->registerCapability(D_COEFFICENT_H2, getComponentName());
kdorheim marked this conversation as resolved.
Show resolved Hide resolved
// Register inputs accepted. Note that more than one component
// can accept an input
core->registerInput(D_EMISSIONS_CO, getComponentName());
core->registerInput(D_EMISSIONS_NMVOC, getComponentName());
core->registerInput(D_EMISSIONS_NOX, getComponentName());
core->registerInput(D_EMISSIONS_H2, getComponentName());
core->registerInput(D_COEFFICENT_H2, getComponentName());


}

Expand Down Expand Up @@ -100,8 +103,8 @@ void OHComponent::setData(const string &varName, const message_data &data) {
H_ASSERT(data.date != Core::undefinedIndex(), "date required");
NMVOC_emissions.set(data.date, data.getUnitval(U_TG_NMVOC));
} else if (varName == D_EMISSIONS_H2) {
H_ASSERT(data.date != Core::undefinedIndex(), "date required");
H2_emissions.set(data.date, data.getUnitval(U_TG_H2));
H_ASSERT(data.date != Core::undefinedIndex(), "date required");
H2_emissions.set(data.date, data.getUnitval(U_TG_H2));
} else if (varName == D_INITIAL_LIFETIME_OH) {
H_ASSERT(data.date == Core::undefinedIndex(), "date not allowed");
TOH0 = data.getUnitval(U_YRS);
Expand All @@ -117,6 +120,9 @@ void OHComponent::setData(const string &varName, const message_data &data) {
} else if (varName == D_COEFFICENT_NOX) {
H_ASSERT(data.date == Core::undefinedIndex(), "date not allowed");
CNOX = data.getUnitval(U_UNDEFINED);
} else if (varName == D_COEFFICENT_H2) {
H_ASSERT(data.date == Core::undefinedIndex(), "date not allowed");
CH2 = data.getUnitval(U_UNDEFINED);
} else {
H_THROW("Unknown variable name while parsing " + getComponentName() +
": " + varName);
Expand Down Expand Up @@ -149,6 +155,8 @@ void OHComponent::run(const double runToDate) {
unitval current_nox = NOX_emissions.get(runToDate);
unitval current_co = CO_emissions.get(runToDate);
unitval current_nmvoc = NMVOC_emissions.get(runToDate);
unitval current_h2 = H2_emissions.get(runToDate);


// get this from CH4 component, this is last year's value
const double previous_ch4 =
Expand All @@ -169,7 +177,13 @@ void OHComponent::run(const double runToDate) {
CNMVOC *
((1.0 * +current_nmvoc) -
NMVOC_emissions.get(NMVOC_emissions.firstdate()).value(U_TG_NMVOC));
toh = a + b + c + d;
const double f =
CH2 *
((1.0 * +current_h2) -
H2_emissions.get(H2_emissions.firstdate()).value(U_TG_H2));


toh = a + b + c + d + f;
H_LOG(logger, Logger::DEBUG)
<< "Year " << runToDate << " toh = " << toh << std::endl;
}
Expand Down Expand Up @@ -201,9 +215,12 @@ unitval OHComponent::getData(const std::string &varName, const double date) {
H_ASSERT(date != Core::undefinedIndex(),
"Date required for NMVOC emissions");
returnval = NMVOC_emissions.get(date);
} else if (varName == D_COEFFICENT_H2) {
H_ASSERT(date == Core::undefinedIndex(), "Date not allowed for H2 coefficent");
returnval = CH2 ;
} else if (varName == D_EMISSIONS_H2) {
H_ASSERT(date != Core::undefinedIndex(), "Date required for H2 emissions");
returnval = H2_emissions.get(date);
H_ASSERT(date != Core::undefinedIndex(), "Date required for H2 emissions");
returnval = H2_emissions.get(date);
} else {
H_THROW("Caller is requesting unknown variable: " + varName);
}
Expand Down
6 changes: 6 additions & 0 deletions src/rcpp_constants.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1503,3 +1503,9 @@ String FLUX_INTERIOR() { return D_FLUX_INTERIOR; }
//' @export
// [[Rcpp::export]]
String HEAT_FLUX() { return D_HEAT_FLUX; }


//' @describeIn oh coefficent for the h2 emissions
//' @export
// [[Rcpp::export]]
String COEFF_H2() { return D_COEFFICENT_H2; }
17 changes: 9 additions & 8 deletions tests/testthat/test-h2.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,9 @@ test_that("tau ", {

test_that("H2 emissions ", {

# Check to make sure that can fetch and set the H2 emissions although
# at this point changing the emissions will have no impact on
# [ch4] dynamics but in the future it should...
# Check to make sure that can fetch and set the H2 emissions and
# that the effects of adding some H2 emissions change thigns
# now that we have the CH4 indrect effects implement.
kdorheim marked this conversation as resolved.
Show resolved Hide resolved

hc <- newcore(ini)
run(hc)
Expand All @@ -45,12 +45,13 @@ test_that("H2 emissions ", {

diff <- abs(out$value - out2$value)

# As of now the [CH4] and tau oh should not change
expect_equal(mean(diff[out$variable == CONCENTRATIONS_CH4()]), 0)
expect_equal(mean(diff[out$variable == LIFETIME_OH()]), 0)
# Now we do expect to see changes in CH4 concentrations in response
# to H2 emissions
expect_gt(mean(diff[out$variable == CONCENTRATIONS_CH4()]), 0)
expect_gt(mean(diff[out$variable == LIFETIME_OH()]), 0)

# But if we can change the H2 emissions we should see a difference in
# H2 emissions between the two new runs!
# Since we are working off of some default H2 emissions = 0, then
# changing the H2 emissions we should see a difference.
expect_equal(mean(diff[out$variable == EMISSIONS_H2()]), new_val)

})
Expand Down
Loading