Update Contaminant Structure Model

[CansuUluseker]

No description provided.

[samharrison7]

I fixed the segmentation fault my removing the line which set the constants variable porosity to 0. It was trying to set this variable before the porosity array was allocated, and in any case, this variable doesn't need a default value. I'm now getting a Failed to read earthworm_densities namelist error that I will try to debug.

[samharrison7]

I've fixed a couple of namelist errors in the latest commit. I'm now getting Given indices do not match output variable rank!. This is coming from the mo_netcdf library, but is being triggered with this stack trace:

#3  0x558f40e81339 in __mo_netcdf_MOD_getreadshape
        at vendor/mo_netcdf/src/mo_netcdf.f90:1869
#4  0x558f40e835b1 in __mo_netcdf_MOD_getdata3df64
        at vendor/mo_netcdf/src/mo_netcdf.f90:1793
#5  0x558f40dd2614 in __datainputmodule_MOD_readbatchvariablesdatabase
        at ././src/Data/DataInputModule.f90:425
#6  0x558f40df4054 in __datainputmodule_MOD_initdatabase
        at ././src/Data/DataInputModule.f90:281
#7  0x558f40d394cd in MAIN__
        at src/main.f90:63
#8  0x558f40d3a7c1 in main
        at src/main.f90:19

I haven't started digging into this yet and have to give up for today, but it looks to be something to do with the emissions data parsing.

[samharrison7]

It looks like you've copied the constants defaults from DefaultsModule to here. Is this a work-in-progress, as I don't think this new file is used anywhere yet? I think it's probably a good idea to separate config and constants defaults, so I'm happy to have these in separate files!

[samharrison7]

This new module is super! 🙂 However, it would be good to add some comments describing each of the procedures and what the code within them is doing 😉

[samharrison7]

Just to check my understanding of this: When you add two contaminants, does this work such that the properties (rho, k_diss etc) stay as the first contaminant in the addition. E.g. for two contaminants A and B, then A + B would result in a contaminant with the properties of A? I think that's fine, but it does make it non-commutative and so we definitely want to document this clearly somewhere! I can imagine my future self being caught out by this 😉

[samharrison7]

This is more a note for the future: A lot of the stuff in this module is specific to particulate contaminants, such as heteroaggregation, so we might want to think about creating a ParticulateContaminant module that extends from this. Other stuff like "pristine", "transformed" and "dissolved" might not be relevant to all contaminants either. I think it's okay to leave this stuff in here for the moment, but this is something to bear in mind when we start adding different contaminants.

[samharrison7]

P.S. This will raise a tonne of difficulties (which is why I suggest not tackling it straight away!). For example, if we go with a ParticulateContaminant module that extends from this module, then we'll probably end up which a rank mismatch for the concentration property c: In this module, c will probably want to be a scalar, whilst in the ParticulateContaminant module, it will at least need a particle size dimension (and then for nanomaterials, it will be the full 3-dimension array that we've got here). There are ways around this, but we need to do a bit of thinking about what would work best.

[samharrison7]

Not sure I understand why this is needed here (or if it is used). Would using the version from the UtilModule cause some kind of dependency loop?

[samharrison7]

Another note for the future: Following on from my comment on ContaminantModule about whether a generic contaminant should have potentially contaminant- (or particulate-) specific processes like attachment, we might want to think about moving process rate constant calculation stuff like this to specific contaminant classes. Again, not something to tackle now, but just bear in mind.

[samharrison7]

(Yet) another option would be to have Reactors specific to different contaminants. In fact, I think perhaps that makes more sense. Either way, it's great to see contaminant-specific code like this ported away from compartments (reaches etc) and into here!

[samharrison7]

Awesome 🙂 Can this procedure be deleted now?

[samharrison7]

Should this file be .gitignored?

[samharrison7]

What's wrong with me?!? 😉 Joking aside, we should probably make this standard across the code. There's a bit of history here: We went with me originally as that was Steve's preference (and I think it comes from VB). I would have at that point in time gone with this instead, mainly because I was writing a lot of Java(Script) and PHP back then. On reflection, I'm wondering whether self is more appropriate for Fortran, based loosely on the fact that Fortran is definitely not C(++), and that Python draws a lot on Fortran roots. Maybe we should put it to vote?! 😆

[samharrison7]

There's a similar argument to be had about whether we use camelCase or snake_case. Again, I went with camelCase originally because I was writing a lot of PHP at the time, but I've come to dislike camelCase (I think it's less readable) and it really doesn't make sense to use something that relies on case sensitivity when Fortran is case insensitive! Obviously changing this would require a full refactor and we probably shouldn't dive into that right now!

[samharrison7]

Nice - I like the backwards compatibility 🤩

[samharrison7]

The latest error is being caused by DATASET%soilBulkDensity not being allocated in NetCDFOutputModule.f90. As far as I can tell, the spatial soil variables (bulk density, water content, texture etc) are no longer in DataInputModule.f90. Have these been moved somewhere or have they got misplaced along the way? I'm pretty sure this is what is causing this error - the model is looking to use variables that haven't been allocated.

This is where they used to be:

nanofase/src/Data/DataInputModule.f90

Lines 433 to 435 in 5a52deb

	! Soil bulk density [kg/m3]
	var = me%nc%getVariable('soil_bulk_density')
	call var%getData(me%soilBulkDensity)

[samharrison7]

This if clause effectively makes quickflow (and runoff below) optional, with a default of 0. I think this is okay, but the docs need updating to reflect this as well (this file: https://github.com/NERC-CEH/nanofase/blob/develop/docs/users/parameters/met-hydro-parameters.csv)

[samharrison7]

(Same for precip)

[samharrison7]

Mostly notes to myself - I made a very brief start on debugging differences in model outputs but will pick this up again after the weekend.

I'm doing it using notebooks (well, a notebook so far) in this repo: https://github.com/samharrison7/fase-output-testing. I've only looked at non-contaminant vars (e.g. SPM concentrations). Findings so far:

• None of the SPM output variables match between the old model (the latest version in the develop branch) and the new model (this PR). This seems to be because, in the new model, SPM concentrations are always 0. Has the input data changed, or is there another issue somewhere?
• Water properties (volume, depth, flow, bed area) do match (within tolorences)
• The dimensions of soil bulk density has flipped between versions. The old model ((y, x)) is the correct way around.

[samharrison7]

I've made a bit of progress debugging. So far I've found that:

• Soil erosion and bank erosion are always 0, when they should be larger than this.
• For soil erosion, this is because the sediment transport capacity is always 0, when it should be larger than this (I haven't looked into bank erosion yet)
• I've only had a quick glance, but I think sediment transport capacity is always 0 because the spatial DATASET%sedimentTransport_a, DATASET%sedimentTransport_b and DATASET%sedimentTransport_c parameters aren't being set in DataInputModule.f90. The model should either take a spatial sediment_transport_a (etc) variable in the NetCDF file, a constant variable in the constant namelist, or default in that in DefaultsModule.f90. However, the spatial variable isn't being checked for now, which means that (when it isn't present) it isn't being defaulted to the constant/default value, and then there must be something setting it to 0 somewhere. These of the relevant lines of code (from the old model) that are missing:

  nanofase/src/Data/DataInputModule.f90

  Lines 545 to 568 in fbc0ae2

  	! Sediment transport param a
  	if (me%nc%hasVariable('sediment_transport_a')) then
  	var = me%nc%getVariable('sediment_transport_a')
  	call var%getData(me%sedimentTransport_a)
  	else
  	allocate(me%sedimentTransport_a(me%gridShape(1), me%gridShape(2)))
  	me%sedimentTransport_a = me%sedimentTransport_aConstant
  	end if
  	! Sediment transport param b
  	if (me%nc%hasVariable('sediment_transport_b')) then
  	var = me%nc%getVariable('sediment_transport_b')
  	call var%getData(me%sedimentTransport_b)
  	else
  	allocate(me%sedimentTransport_b(me%gridShape(1), me%gridShape(2)))
  	me%sedimentTransport_b = me%sedimentTransport_bConstant
  	end if
  	! Sediment transport param c
  	if (me%nc%hasVariable('sediment_transport_c')) then
  	var = me%nc%getVariable('sediment_transport_c')
  	call var%getData(me%sedimentTransport_c)
  	else
  	allocate(me%sedimentTransport_c(me%gridShape(1), me%gridShape(2)))
  	me%sedimentTransport_c = me%sedimentTransport_cConstant
  	end if

  .

I'll try and dig a bit deeper tomorrow!

[samharrison7]

I'm making progress! I re-added the sediment calibration parameters (which had disappeared from DataInputModule.f90): deposition_alpha, deposition_beta, resuspension_alpha, resuspension_beta, sediment_transport_a, sediment_transport_b, sediment_transport_c, bank_erosion_alpha and bank_erosion_beta

Now the SPM concs are close but the new predictions are always higher than the old. Here is the SPM conc timeseries for all of the grid cells to show what I mean:

[samharrison7]

Think it's something to do with resuspension as this is always zero in the updated model. The resus params are being passed in correctly, but something must be going wrong with the resus rate calculation.

[samharrison7]

Latest progress: Resuspension parameters are being read in correctly the resuspension rates (k_resus) seem to being set correctly. The bit that is causing resuspension to be zero is that Mf_bed_by_size() returns odd results (used in updateDisplacementRiverReach).

In the new model, Mf_bed_by_size() is constant and mostly zero:

mf_bed_by_size   0.0000000000000000        0.0000000000000000        81.999998167157187        0.0000000000000000        0.0000000000000000

In the old model, the first size class varies, e.g.:

 mf_bed_by_size   2.4828609780654038E-003   8.0000000000000002E-002   1.2800000000000000        1.2800000000000000        1.2800000000000000
 mf_bed_by_size   2.4828824263178630E-003   8.0000000000000002E-002   1.2800000000000000        1.2800000000000000        1.2800000000000000
 mf_bed_by_size   2.4828794795080937E-003   8.0000000000000002E-002   1.2800000000000000        1.2800000000000000        1.2800000000000000
 mf_bed_by_size   5.4973661057711462E-003   8.0000000000000002E-002   1.2800000000000000        1.2800000000000000        1.2800000000000000
 mf_bed_by_size   5.5720733442795054E-003   8.0000000000000002E-002   1.2800000000000000        1.2800000000000000        1.2800000000000000

I have to leave this now, but I've pushed my changes in the latest commit.

[samharrison7]

Almost a year later(!). I was still struggling with sediment issues, but after some debugging have hopefully managed to fix it. The culprit was that distributionSediment was only being derived if no soil data was present, as it had moved inside the else block of a new if (haveSoil2D) check. This meant that, if 2D soil data are present, distributionSediment defaulted to the default_spm_size_distribution (50/30/10/7/3) placeholder and soil texture wasn't used to calculate the distribution. Applying this fix means SPM concs now match (see below). I've also applied this fix on the PFAS branch, which I'll push soon.

[samharrison7]

I'll hold off merging this as it might make more sense to merge in the PFAS modifications, which (I think?) contain all the contaminant-agnostic changes here.