Changed a solver for a soil moisture reservoir to prevent over-extraction of water by primary fluxes

[RY4GIT]

[Short description explaining the high-level reason for the pull request]

Current behavior

• There was an issue in the soil conceptual reservoir: the secondary fluxes was constantly zero. This was because the primary flux was over-extracting water from the soil reservoir, and there was no water left for the secondary fluxes.
• The root cause of the problem was in the numerical solver: the model before this commit was extracting fluxes from the soil reservoir in order (ET -> primary flux -> secondary flux), and the secondary fluxes are calculated based on soil water storage remaining after ET and primary flux used up the soil water storage. This "extracting in order" method does not work unless the time step is super fina e (like 1-minute interval or so).

Expected behavior

• The solution is to use ODE for the soil reservoir. With this, the fluxes from the soil reservoir (ET, primary, and secondary flux) are extracted simultaneously with a finer timestep (but not too fine). Adding this ODE slows the calculation but fixes the unrealistic soil moisture reservoir behavior.

Additions

• Added an ODE solving module for the soil moisture conceptual reservoir calculation. The main functions/methods added are: conceptual_reservoir_flux_calc, jac, and soil_reservoir_flux_calc.

Changes

• Some variables are added or changed according to this change (e.g. diff_infilt, diff_perc, )

Deletions

• A few lines of code that compares groundwater water deficit and percolation are removed (the ones starting from cfe_state.previous_flux_perc_m > cfe_state.soil_reservoir_storage_deficit_m), because the ODE will solve the issue of excess percolation without those lines
• A method et_from_soil is removed, and incorporated into the soil_reservoir_flux_calc

Testing

1. Visually checked plots, checked mass balance by running run_massbalance_cfe.py (water balance is fixed by my commit on Sept 22th)

Notes

• Calculation time, unfortunately, slows down according to this change. My test case took 30 seconds to calculate 3 years of record. To minimize the calculation time, I tested different solvers and picked the best one, as well as added the jacobian matrix.
• Currently the code uses the scipy odeint for ODE solver, which uses the LSODA algorithm. The rationale behind it is that (1) an ODE solver for stiff system was needed as the target system (soil moisture reservoir) can be unstable around the thresholds (2) an ODE solver that does not take too much computational time was preferred; LSODA switches between stiff and non-stiff methods, which saves computational time. Note that there are potential ODE solver options, such as to use solve_ivp, and try different ODE solver scheme for better accuracy in water balance.

Checklist

• PR has an informative and human-readable title
• Changes are limited to a single goal (no scope creep)
• Code can be automatically merged (no conflicts)
• Code follows project standards (link if applicable)
• Passes all existing automated tests
• Any change in functionality is tested
• New functions are documented (with a description, list of inputs, and expected output)
• Placeholder code is flagged / future todos are captured in comments
• Visually tested in supported browsers and devices (see checklist below 👇)
• Project documentation has been updated (including the "Unreleased" section of the CHANGELOG)
• Reviewers requested with the Reviewers tool ➡️

Target Environment support

• Windows
• Linux
• Browser