Add new feature: Lateral flow of Dissolved Organic Carbon (DOC)

[devarajun]

To begin with, solving the mass balance equation for DOC transport is planned. This will imply changes in MOSART (https://github.com/ESCOMP/MOSART/blob/master/src/riverroute/MOSART_physics_mod.F90)
and in CTSM Soil Carbon Module (https://github.com/ESCOMP/CTSM/blob/master/src/biogeochem/CNVegCarbonStateType.F90).

Brief science information:

Assume the lateral transport of DOC to be proportional to the water fluxes, i.e. the export from one reservoir to the next have the same concentration of DOC. DOC is drawn up into runoff from the top 5cm of soil layer.

At each spatial unit (lat/lon grid), DOC goes into water fluxes generated from hillslopes, then tributaries and to main channel. DOC implementation code mainly builds upon on MOSART code for taking up DOC into hillslope, tributaries and main channel.

Plan for the implementation of the DOC mass balance code:

We plan to implement the changes for Soil Carbon module and MOSART separately in their respective repositories.

Pseudocode for MOSART

DOC mass balance equation mainly requires the total water storage for each category of hydro-logical unit and the channel flow rate. This information used in the mass balance equation of DOC transport.

Pseudocode for CTSM Soil carbon module

DOC production in soils derived from the first-order kinetics formulation for different carbon pools (see the article for the DOC production formula https://gmd.copernicus.org/articles/11/593/2018/).

Long term goal

• Long term plan is to couple the microbial soil decomposition model (developed by Elin Ristorp Aas, PhD Student at University of Oslo) with MOSART-DOC.

• If possible to extend simple DOC mass balance to follow Liao et al (2019) https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001792

People involved: @devarajun , @sunnivin , @ecaas (Elin Ristorp Aas) @kjetilaas (Kjetil Schanke Aas)

[wwieder]

Thanks for filing this issue, @devarajun. Looking forward to seeing your contributions down the road.
It may be worth coordinating the efforts to integrate the land changes to soil BGC with mizuRoute #920 instead (or in addition) to MOSART? Also, I'm assuming you'll produce DON fluxes with the DOC fluxes? We may as pass these to the river model too?

[devarajun]

Thank you @wwieder for commenting. Regarding DON, yes we have plans for DON fluxes. We were not aware of mizuRoute and thanks for pointing this out. How different mizuRoute from MOSART in terms of biases in high latitudes? I will try to explore mizuRoute documentation. Is it already a optional component of CTSM/CLM and works in regional CLM simulation setup at any resolution?

[wwieder]

I think don't think mizuRoute is ready to go yet, but an upcoming feature that's being actively developed now. @nmizukami and @ekluzek can you confirm? @devarajun I mainly wanted to give you a heads up that mizuRoute is coming. I trust adding tracer fluxes aren't an overwhelming burden in either river model.

[dlawrenncar]

Minor point, but should stress that we don't really expect streamflow biases to be strongly affected by using MOSART or mizuRoute, especially for big rivers. The primary source of the biases in most cases, I think, is the amount of runoff that is generated by CTSM. There are other benefits to mizuRoute beyond reducing biases though, including better integration between rivers and lakes and a more accurate river network, especially for smaller rivers.

…

On Thu, Nov 19, 2020 at 5:29 AM will wieder ***@***.***> wrote: I think don't think mizuRoute is ready to go yet, but an upcoming feature that's being actively developed now. @nmizukami <https://github.com/nmizukami> and @ekluzek <https://github.com/ekluzek> can you confirm? @devarajun <https://github.com/devarajun> I mainly wanted to give you a heads up that mizuRoute is coming. I trust adding tracer fluxes aren't an overwhelming burden in either river model. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1216 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFABYVG4BAGDNWSSVVBLJRLSQUFUHANCNFSM4TYZ2Z3A> .

[ekluzek]

@wwieder we do have a preliminary version of mizuRoute that comes with CTSM as of ctsm5.1.dev016. It doesn't really function there though, but you can at least look at a version of the code. I have a branch on my fork add_mizuRoute that we have working for a few configurations. We'll likely only have mizuRoute working with CTSM for a limited set of resolutions.

Also mizuRoute isn't currently handling tracers. @nmizukami, and @martynpclark will have to let us know about how hard that will be to add. As @dlawrenncar points out it might make the most sense to add certain features to only one river model (especially at first) for the configuration that it makes the most sense in. Right now it seems like the roles are: RTM will be the river model for Paleo work, MOSART for climate, and mizuRoute for NWP. Dynamic lakes is a feature that's just going into mizuRoute, so we already have precedence for features going into one and not all river models. And practically you can't always put all needed features into each of the different river models. I imagine in the long term it's possible mizuRoute will replace both RTM and MOSART, but right now we need all three.

[nmizukami]

Hi @devarajun, adding to comments from @ekluzek, the main difference between mizuroute and mosart/rtm is river network. mosart/rtm use gridded network while mizuRoute uses catchments (you can think this as unstructured mesh, but complex shape). The difference in river network can change river discharge if drainage area is different.
For tracer capability, I think we (including @martynpclark) will need to discuss the best path forward.

[ekluzek]

Another aspect of this is that DOC will need to be passed from CTSM through the coupler to MOSART. We'll likely want to implement this in the NUOPC coupler CMEPS rather than the current MCT coupler. There might also need to be an XML variable to trigger that this should be turned on in both CTSM and MOSART. This would likely be similar to this issue for the support of mizuRoute...

ESCOMP/CMEPS#132

[mvertens]

@ekluzek - thanks for pointing out the implementation in CMEPS. I am happy to help with this and review the implementation.

[devarajun]

Thanks to all of you for your useful inputs. @ekluzek great that you have pointed out on passing DOC through the coupler and I now know this at the beginning of coding thanks! @nmizukami @mvertens will keep in touch with you.

[devarajun]

Hi, @ekluzek @nmizukami @mvertens @wwieder @dlawrenncar
As I explained earlier my intention was to implement a subroutine directly into MOSART. However, considering mizuRoute (and RTM) river routing options you have mentioned, if its feasible to make DOC/DON work with any current and future river routing. What do you think is the best way of implementing DOC/DON code? It would be straight forward to include a DOC/DON subroutine in either of the river routing components. But with respect to future developments and compatibility, I wonder if creating a module or subroutine for DOC/DON in CLM biogeochem source code that access stream-flow in & out and runoff through the coupler from either MOSART or mizuRoute, could be a better choice? This would require or at least benefit from a well defined interface to river routing sub-components. By the way, the mass-balance equation of tracers along with the source term given below mainly requires lateral flow in and out for each hydrological unit and runoff generation (for eg. hillslope, sub-network, and the main channel for MOSART). Q→ Streamflow, c→ concentration of DOC/DON, r→ DOC source/sink term.
d( CS)/dt = sum(QC)_inflows - sum (Q*C)_outflows + r

[nmizukami]

It seems to me that it would not be very difficult to implement the tracer based on mass balance equation like that (advection only) into river models (without any details in each river model).

I wonder if creating a module or subroutine for DOC/DON in CLM biogeochem source code that access stream-flow in & out and runoff through the coupler from either MOSART or mizuRoute, could be a better choice?

Your routine needs to know flow direction to move DOC/DON downstream at land model grid? IMO: DOC/DON is passed from land model to river, and the tracer subroutine should reside in the river model component.

Which river models should be used? maybe others can weigh in?

[dlawrenncar]

For now, I think the main river model should be MOSART with the goal of also eventually being able to use mizuRoute. I don't know enough to be able to anticipate whether or not it would be possible to make this general enough to be applicable across multiple river models. That will take some thought and exploration. Perhaps we should meet after the holidays to discuss.

…

On Tue, Dec 22, 2020 at 12:08 PM Naoki Mizukami ***@***.***> wrote: It seems to me that it would not be very difficult to implement the tracer based on mass balance equation like that (advection only) into river models (without any details in each river model). I wonder if creating a module or subroutine for DOC/DON in CLM biogeochem source code that access stream-flow in & out and runoff through the coupler from either MOSART or mizuRoute, could be a better choice? Your routine needs to know flow direction to move DOC/DON downstream at land model grid? IMO: DOC/DON is passed from land model to river, and the tracer subroutine should reside in the river model component. Which river models should be used? maybe others can weigh in? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1216 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFABYVBKP6CAK6H7N37C7KTSWDVEPANCNFSM4TYZ2Z3A> .

[devarajun]

Thanks @nmizukami @dlawrenncar for the quick response. Yes my understanding is also that DOC/DON should reside in the river model component in order to route tracer along with discharge. I am happy to meet and discuss sometimes in January if there is any possibility to make this general enough to be applicable to different river models.

[wwieder]

I don't have any other suggestions, other than the naive assumption that once you pass the C/N fluxes to the river it's just building the coupler infrastructure to pass them to the river model and as @ekluzek mentioned this should be done in NUOPC. Also agree with @dlawrenncar, let's have a call in 2021 to discuss.

[mariuslam]

As @devarajun has filled a new position, I am now supposed to take over.

Thank you @wwieder, @nmizukami , @ekluzek, for the messages above and in issue (#1458), they have been really helpful to get DOM fluxes into the existing structure of the decomposition cascade, by tracking down the way it is done for heterotrophic respiration (HR).

• We have successfully passed DOC from CTSM to MOSART through CMEPS (branches are listed below), by tracking the way it is done for Flrl_rofsur.
• In MOSART, we added a dimension to nt_rtm: "nliq", "nice" and now also "ndom". We use ndom to send DOM to the variable "rtmCTL%qsur()", which is the input surface runnof.
• We added history fields, such as the discharge_to_ocean and storage (shown below).

But after discussing with @devarajun, there are some things we did not figure out yet.

1. We assume that the outputted DOM variables in MOSART currently have the wrong units. DOM is implemented as a tracer the same way liquid water is but the units of DOM in CTSM are gC/m2/s instead of mmH20/s for water.
2. We want DOM to be transported/distributed proportionally to the amount of liquid water, is this currently the case ? Or should there be some additional mass balance equations for DOM?

Thank you @wwieder, @nmizukami and others for suggestions, any inputs are welcome. Also I haven't looked at Mizuroute but may do so if recommended.

Further down the road we are planning to better parametrize the production rates of DOM, and validate at some sites. I am eager to collaborate on these tasks with whoever has some input :)

• https://github.com/MetOs-UiO/CTSM/tree/dom
• https://github.com/MetOs-UiO/CMEPS/tree/dom_tracer
• https://github.com/MetOs-UiO/MOSART/tree/dom_tracer

[wwieder]

thanks for renewing this conversation @mariuslam. Glad you're going to continue this work. What do you need from us to keep making progress? I have to confess, once DOM leaves land I don't know what to think about it, although it makes sense to convert units to more meaningful aquatic units (fluxes = gC/L/s, concentrations = gC/L?)

One notable shift is that we'll need to make sure the code works with the nuopc coupler, not mct (which was used when the project started). Towards this end, I wonder how hard it is to bring your development branch(s) up to a more modern CTSM5.1 development branch?

On a more technical note, can we specify if we're talking about DOC or DON (instead of DOM)?

[nmizukami]

Hi @mariuslam, thanks for posting this.

..... DOM is implemented as a tracer the same way liquid water is....

2. We want DOM to be transported/distributed proportionally to the amount of liquid water, is this currently the case ? Or should there be some additional mass balance equations for DOM?

Looking through your changes in MOSART, i don find where DOM is transported through network using flow velocity. MOSART has two tracers (liquid and ice) and you added dom. liquid is moved using Euler subroutine in MOSART_physics_mod.F90 (it uses kinematic wave). For Ice, i believe it is moved to directly outlet here (correct me if this is wrong). I don't see how dom tracer is handled after it is imported into MOSART. So wondering how TOTAL_DISCHARGE_TO_OCEAN_DOM is computed. i might miss something.

If dom is not moving through network with flow velocity, this will have to be done I believe. I was thinking of using advection-diffusion equation (equivalent to diffusive wave equation in river routing, except variable is concentration instead of discharge) or starting with only advection, which is equivalent to kinematic wave equation used in MOSART euler. so you could move dom downstream with flow velocity computed in MOSART (was this done??)

[mariuslam]

Hi, Thank you for the quick answers.
@wwieder, good point, I will rename the variable to DOC, so we can add DON later on. Concerning the CTSM branch, it is relatively recent (ctsm5.1.dev112), and I use the NUOPC coupler.
@nmizukami, thank you for having a quick look at the changes! My understanding is that DOM becomes nt_rtm=3, so unlike for ice, Euler_calc remains equal to TRUE
(https://github.com/MetOs-UiO/MOSART/blob/064a31de847460e49430fa51a27388b7ac3543b9/src/riverroute/MOSART_physics_mod.F90#L59)
So I think that there is no need to add code if we want DOM to be transported with flow velocity the same way water is, (using only advection then i guess). Does it make sense or is this is wrong? Besides the inaccurate units of the outputted variables, I wonder if we can just forget about converting DOM to concentrations during the transport (solving of Euler equations)? Does flow velocity results in DOM being transported proportionally to the water volume?

I don't know if any other MOSART user, not yet in the conversation may have some input ? (I see Hongyi Li has made some functions of mosart).

Thanks for asking @wwieder, I posted the message to give a little update, see if there was no obvious mistake yet, but I hope I can still make some progress by exploring a bit on my own for the moment.

[nmizukami]

Hi @mariuslam, thanks for pointing me to L59 in MOSART_physics_mod.F90. Yes, DOM tracer (nt_rtm=3) does go through Euler subroutine. However, I still doubt that using Euler subroutine directly to move DOM as is is right.

If you look through all the subroutines used in the Euler routine, the Euler routine computes flow velocity using the Manning equation (function CRVRMAN). This is applicable to only water (it uses a hydraulic variable- hydraulic radius and roughness coeff. of the river bed, and hydraulic radius is computed using water storage, and flow area from the previous time step I believe), and I don't think it makes sense to apply the Manning equation to DOM concentration (value computed from the manning eq. with DOM input is physically meaningless). I think answering to your question is no

Does flow velocity results in DOM being transported proportionally to the water volume?

I think you want to use flow velocity computed from liquid tracer for DOM advection.

[mariuslam]

Thanks @nmizukami.
Then @devarajun was right to create 3 additional equations to solve DOM with the Euler algorithm.
My question then is: should I (1) turn Euler_calc to false for DOM and go into my own subroutines but still keep it as a tracer (i.e. meaning nt_rtm=3) and use the same variable names as for LIQ and ICE along the routing network?
Or (2) add new variables to track the flow of DOM, but not add a column to to the existing variables, and keep nt=2 (tracer number)
Cheers, Marius

[nmizukami]

Hi @mariuslam, This is my thought, but need to discuss with others.
For 1), I think it would be cleaner to create a separate routine for a tracer of any constituents, and turn Euler_calc to false so not going to Euler.
For 2), I think maybe need to look at how many variables in Liq and Ice TRunoff are shared with DOM (I don't think many, maybe none?). My guess is use a separate data structure than TRunoff data structure.

[mariuslam]

Ok,It will be cleaner with a new data structure for DOM, so I can give more adapted names and units to the variables.
But I guess we want to avoid having TRunoff variables with wrong/empty values in the third column (nt_rtm=3, DOC), and later on maybe in 4th collumn DON (nitrogen).
Would it make sense to keep nt_rtm=2 and instead add DOC into a new tracer variable (nt_dom for example)?
Then if I create a new data structure as you suggest (TDom for example), the variables of Tdom will only have the dimension of nt_dom, so they wont have empty/wrong values in the first two columns as they would have had if I used nt_rtm=3.

I hope you ll understand what I mean :)

[nmizukami]

I think we are on the same page!

[mariuslam]

Hi, @devarajun, @nmizukami, @ekluzek
I have updated the code according to what we discussed earlier. (See: https://github.com/MetOs-UiO/MOSART/tree/dom_tracer_2)

Although the model succeeds, given the output, I am pretty sure there are several issues with the code.
I am for example confused with the units: before Euler is called, water inputs are converted from m3/s to m/s here. Does that mean that domH, domT and domR are in kg/m?

Does this also means that wh, wt, wr are in m here? But why is only wh multiplied by the area then?

I guess there may be other issues than the units, does anyone has tips on how to check that the DOC flow is correct? Right now I am running global as it is the most straightforward.

Thanks, Marius

[mariuslam]

Hi,
The structure for DOM transport by advection in Mosart with DOC as a tracer is now implemented.
I will look into how to better represent the production of DOC in CTSM, for now I use fixed fractions like for respiration, but maybe this can be elaborated? Any other ideas of future improvements before I evaluate the model are welcome @nmizukami , @wwieder :)
Here some preliminary results.

[nmizukami]

Hi @mariuslam,

I guess I missed your earlier comments about units. Not sure if you figured out? MOSART uses a few separate data structure - rtmCTL used for communication between coupler and mosart, and TRunoff for actual mosart physics (used in euler). and all the units are in RunoffMod.F90.

I don't have good sense on the reasonable magnitude of concentration.

I guess it may be easier to evaluate the validity of model outputs by focus some basin and look at time series at some river point (e.g., outlet) and DOC follows discharge? and look at budget?

[mariuslam]

Hi @nmizukami,

Yes thank you, I managed to figure out where units are changed during the passing from the coupler, so the units should be correct now. (If DOC is wrong, it is due to the production rate I assigned in CTSM.)

I have access to several DOC concentration data sources, including some at Norwegian river outlets which will be useful when evaluating. But I may want to work a bit more on the code first.

Here are ideas:

• improve the DOC production parameters in CTSM decomposition cascade?
• what to do with recalcitrant and labile pools ? (now just summed and sent to MOSART as surface and subsurface DOC)
• add DON transport ?
• decomposition of DOM in the river?
• flooding of DOM back to land?
• Add some DOM sink in Ocean model?
• Complexify the transport in Mosart? @nmizukami you suggest advection-diffusion, but I guess this means changing water flow in Mosart and is not directly linked to the DOM mass balance equations already implemented right?

Thanks in advance for feedback on this !

[wwieder]

This is awesome to see. Are you able to present at the LMWG this Feb, @mariuslam ? Registration closes today.
I'm also curious to know what you feel like your timeline is for this work? I wonder if we should support the river transport in MOSART long-term, or really just focus on MIZUROUTE? These larger CESM integration questions may take more time to iron out.

[mariuslam]

Thanks @wwieder. I signed up at LMWG but without talk as I only have preliminary results (no-spinup, no observations). If you judge this should be added to the plan, I could give a 5 min talk.
I will be highly focused on DOM-MOSART until July. I think my main goal will be to evaluate the DOC production and river discharge against observations and adjust the DOC production parameters accordingly (or make this more dynamical). This also seems safer as it is on the CTSM side...
When the CESM integration questions will be ironed out, depending on the timing, I can try to have a look at mizuroute, I don't know how hard it would be to transfer this work there @nmizukami ? Otherwise, a PhD in Oslo may continue developing this further after I leave.

[wwieder]

OK, this may warrant a slightly longer scientific discussion at a Thursday meeting? We'll see how the schedule looks for the LMWG and I'll let you know, @mariuslam . My guess is it would be a valuable conversation with at the BGCWG, as the export of DOC and DON to the ocean is something they may be keen to discuss?

I'm adding next so we can briefly discuss at CTSM SE meeting, but I think the scientific implications of this work should be evaluated first, before deciding on software priorities.

[mariuslam]

Hi,

I have been struggling a bit in the past weeks, trying to figure out why the DOC entering Mosart is 200 times smaller than what is going out to the ocean.

There must be an error in the way I coded this, but I can't find it. I would really appreciate som help or a review of what I did if this is at all possible. https://github.com/MetOs-UiO/MOSART/tree/Dev_dom
@ekluzek , @billsacks maybe if you have time ?

Marius

[wwieder]

Hi @mariuslam, we'll discuss at our SE meeting this week, although at this time I'm afraid we're all stretched really thin to commit much time to this issue.

[nmizukami]

Hi @mariuslam, sorry to hear that and this is frustrating.

maybe it would be easier to see the code changes and review them if you pull request (as a draft). I assume you eventually want to merge into ESCOMP/MOSART? maybe @billsacks and @ekluzek can chime in?

you might have checked this, that sounds like there is something wrong in unit-conversions somewhere?

Is this based on budget check like water here?

[ekluzek]

I do endorse a PR as a nice way to be able to share code with others (so agree with @nmizukami above). We do sometimes open PR's just for visibility even if we know something isn't coming in -- just as a way to make it easy for people to view it.

I can maybe help with some coding issues, but for this we really need someone familiar with river hydrology. So @nmizukami is one of the best resources here, but maybe also @swensosc or Andy Wood? I also agree with @billsacks that I don't have time to devote to this right now either.

[mariuslam]

Hi,
Thanks you all for your answers.
I made a PR as suggested: ESCOMP/MOSART#64
@nmizukami, Yes I have been reading my changes multiple times and checking the units. I added a few if statements because the water content seems to become negative in some grid cells. Maybe it has something to do with that?
It may also be because of the way I extract these variables from Mosart. It is weird that there is a loss of DOC in the hill but a large gain in the subnetwork.
@ekluzek, thank you for the suggestions. In the meanwhile I will continue looking for possible causes.
These are DOC numbers integrated over a year:
585953096.1560751 qsur
2272295019.127152 qsub
499076800.08670616 out_hills
539132943500.8501 out_subn
649837252643.7626 river discharge ocean
These are for water:
11730926098465.37 qsur
3810542169065.1494 qsub
12694138488383.57 out_hills
12591109529608.68 out_subn
15627459302912.102 river discharge ocean

[nmizukami]

Hi @mariuslam,

I think this is good information. so I understood qsur~= out_hills, and qsur+ qsub ~= out_subn = river discharge ocean in MOSART.

I looked at three subroutines in DommasbMod, and I am wondering the parts using min and max functions is causing off-balance. I put the comments on the pull-request.

if you get negative flow (this is probably separate issue), that will cause negative concentrations, but if you restrict domHout >=0 (change concentration), it cause off balance?

what would happen if you comment out max/min functions there and run? you might get negative flux, but does it balance?

[mariuslam]

Hi @nmizukami ,

I think you understand correctly concerning the fluxes that should be equal for correct balance.
I ran the model commenting out all the min and max functions in Dommas.F90 file.

Weirdly I get exactly the same results.
So the problem should be elsewhere in the code, I think this huge jump in subnetwork is weird.

Regards,

Marius