SIP256 Leaf on/off events

.github/workflows/link-check.yml

@@ -30,6 +30,7 @@ jobs:
             --format markdown
             --exclude '^https://doi\.org/'
             --exclude '^https://russellranch\.ucdavis\.edu'
+            --exclude '^https://join\.slack\.com'
             .
           fail: true
         env:

docs/CHANGELOG.md

@@ -27,6 +27,7 @@ sections to include in release notes:
 ### Added
 
 - `sipnet-view` tool for visualizing SIPNET output files (#317)
+- `leafon` and `leafoff` events for tracking phenological transitions (#326)
 
 ### Fixed
 

docs/README.md

@@ -24,7 +24,7 @@ SIPNET (Simplified Photosynthesis and Evapotranspiration Model) is a lightweight
 - **Understanding the code:** The [developer docs](model-structure.md) describe architecture, while sections on [testing](developer-guide/testing.md) and [CLI extensions](developer-guide/cli-options.md) cover contributions.
 - **Project practices:** Review the [Contributing guide](CONTRIBUTING.md), [Code of Conduct](CODE_OF_CONDUCT.md), and [Changelog](CHANGELOG.md) before opening PRs or issues.
 
-Need help? Open an issue on [GitHub](https://github.com/PecanProject/sipnet/issues) or join the [PEcAn community Slack](https://join.slack.com/t/pecanproject/shared_invite/zt-3ile31ylu-r1sMh~esl~7A_AolYNNYLQ).
+Need help? Open an issue on [GitHub](https://github.com/PecanProject/sipnet/issues) or join the [PEcAn community Slack](https://join.slack.com/t/pecanproject/shared_invite/zt-3wlzwms7n-mmoLuYsIfaDqp_FRZD1eQQ).
 
 ## License
 

docs/model-structure.md

@@ -1162,6 +1162,32 @@ f_{\text{intercept}} \, F^W_{\text{irrig}}, & I_{\text{irrigation}} = 0 \\
 \label{eq:irrig_evap}
 \end{equation}
 
+### Leaf On/Leaf Off
+
+Leaf on and leaf off events define the timing of leaf emergence and senescence, respectively. These events directly
+specify the amount of carbon added to the leaf carbon pool on the leaf on date, and the fraction of carbon removed from
+the leaf carbon pool on the leaf off date. 
+
+When a leaf on event occurs, an amount of carbon (specified by the `leafGrowth` parameter) is transferred from the wood
+carbon pool to the leaf carbon pool. As leaf C:N is usually lower than wood C:N, the excess nitrogen
+implied by the static C:N ratios is included as part of the plant nitrogen demand. If there is insufficient nitrogen
+available for this lump-sum move, nitrogen limitation will occur. 
+
+When a leaf off event occurs, a fraction of the leaf carbon (specified by the `fracLeafFall` parameter) is transferred
+from the leaf carbon pool to the litter pool (or soil pool, if the litter pool is not being used). The corresponding
+nitrogen (calculated from the leaf C:N ratio) is also transferred to the litter or soil nitrogen pool.
+
+**Event parameters:**
+
+| Parameter | Value                | Description       |
+|-----------|----------------------|-------------------|
+| Year      | integer              | Year              |
+| Day       | integer              | Day of year       |
+| Type      | `leafon` / `leafoff` | The type of event |
+
+There are no other parameters needed for these events, as the amount of transfer is determined by the parameters
+mentioned above.
+
 <!-- 
 **Flooding** increases soil water to water holding capacity and then adds water equivalent to the depth of flooding. Subsequent irrigation events maintain flooding by topping off water content.
 

docs/parameters.md

@@ -154,21 +154,21 @@ Run-time parameters can change from one run to the next, or when the model is st
 
 ### Phenology-Related Parameters
 
-| Symbol                   | Parameter Name   | Definition                                                            | Units                            | Notes                                |
-| ------------------------ | ---------------- | --------------------------------------------------------------------- | -------------------------------- | ------------------------------------ |
-| $D_{\text{on}}$          | leafOnDay        | Day of year when leaves appear                                        | unitless                         | day of year (1–365)                  |
-| $GDD_{\text{on}}$        | gddLeafOn        | GDD threshold for leaf appearance (GDD-based phenology)               | $°\text{C} \cdot \text{day}$     |                                      |
-| $T_{\text{on}}$          | soilTempLeafOn   | Soil temperature threshold for leaf appearance (temp-based phenology) | $°\text{C}$                      |                                      |
-| $D_{\text{off}}$         | leafOffDay       | Day of year for leaf drop                                             | unitless                         | day of year (1–365)                  |
-| $\Delta C_{\text{leaf}}$ | leafGrowth       | Additional leaf growth at start of growing season                     | $\text{g C} \cdot \text{m}^{-2}$ |                                      |
-| $f_{\text{fall}}$        | fracLeafFall     | Additional fraction of leaves that fall at end of growing season      | unitless                         |                                      |
-| $\alpha_{\text{leaf}}$   | leafAllocation   | Fraction of $NPP$ allocated to leaf growth                            | unitless                         |                                      |
-| $K_{\text{leaf}}$        | leafTurnoverRate | Average turnover rate of leaves                                       | $\text{year}^{-1}$               | Converted to per-day rate internally |
+| Symbol                    | Parameter Name   | Definition                                                            | Units                            | Notes                                |
+|---------------------------|------------------|-----------------------------------------------------------------------|----------------------------------|--------------------------------------|
+| $D_{\text{on}}$           | leafOnDay        | Day of year when leaves appear                                        | unitless                         | day of year (1–365); 0 to turn off   |
+| $GDD_{\text{on}}$         | gddLeafOn        | GDD threshold for leaf appearance (GDD-based phenology)               | $°\text{C} \cdot \text{day}$     |                                      |
+| $T_{\text{on}}$           | soilTempLeafOn   | Soil temperature threshold for leaf appearance (temp-based phenology) | $°\text{C}$                      |                                      |
+| $D_{\text{off}}$          | leafOffDay       | Day of year for leaf drop                                             | unitless                         | day of year (1–365); 0 to turn off   |
+| $\Delta C_{\text{leaf}}$  | leafGrowth       | Additional leaf growth at start of growing season                     | $\text{g C} \cdot \text{m}^{-2}$ |                                      |
+| $f_{\text{fall}}$         | fracLeafFall     | Additional fraction of leaves that fall at end of growing season      | unitless                         |                                      |
+| $\alpha_{\text{leaf}}$    | leafAllocation   | Fraction of $NPP$ allocated to leaf growth                            | unitless                         |                                      |
+| $K_{\text{leaf}}$         | leafTurnoverRate | Average turnover rate of leaves                                       | $\text{year}^{-1}$               | Converted to per-day rate internally |
 
 ### Allocation Parameters
 
 | Symbol                        | Parameter Name       | Definition                                  | Units    | Notes                                                                                                                                             |
-| ----------------------------- | -------------------- | ------------------------------------------- | -------- | ------------------------------------------------------------------------------------------------------------------------------------------------- |
+|-------------------------------|----------------------|---------------------------------------------|----------|---------------------------------------------------------------------------------------------------------------------------------------------------|
 | $\alpha_{\text{fine root}}$   | fineRootAllocation   | Fraction of $NPP$ allocated to fine roots   | unitless |                                                                                                                                                   |
 | $\alpha_{\text{coarse root}}$ | coarseRootAllocation | Fraction of $NPP$ allocated to coarse roots | unitless | Calculated internally from remainder: $\alpha_{\text{coarse root}} = 1 - \alpha_{\text{leaf}} - \alpha_{\text{wood}} - \alpha_{\text{fine root}}$ |
 | $\alpha_{\text{wood}}$        | woodAllocation       | Fraction of $NPP$ allocated to wood         | unitless |                                                                                                                                                   |

src/sipnet/events.c

@@ -136,6 +136,32 @@ EventNode *createEventNode(int year, int day, int eventType,
       tParams->tillageEffect = tillEffect;
       newEvent->eventParams = tParams;
     } break;
+    case LEAFON: {
+      double dummy;
+      LeafOnParams *lParams = (LeafOnParams *)malloc(sizeof(LeafOnParams));
+      // Check for extraneous data: leafon takes no parameters, so error if any
+      // numbers are found. sscanf returns 0 or EOF (-1) for empty/whitespace.
+      int numRead = sscanf(eventParamsStr,  // NOLINT
+                           "%lf", &dummy);
+      if (numRead > NUM_LEAFON_PARAMS) {
+        logError("parsing LeafOn params for year %d day %d\n", year, day);
+        exit(EXIT_CODE_INPUT_FILE_ERROR);
+      }
+      newEvent->eventParams = lParams;
+    } break;
+    case LEAFOFF: {
+      double dummy;
+      LeafOffParams *lParams = (LeafOffParams *)malloc(sizeof(LeafOffParams));
+      // Check for extraneous data: leafoff takes no parameters, so error if
+      // any numbers are found. sscanf returns 0 or EOF (-1) for empty input.
+      int numRead = sscanf(eventParamsStr,  // NOLINT
+                           "%lf", &dummy);
+      if (numRead > NUM_LEAFOFF_PARAMS) {
+        logError("parsing LeafOff params for year %d day %d\n", year, day);
+        exit(EXIT_CODE_INPUT_FILE_ERROR);
+      }
+      newEvent->eventParams = lParams;
+    } break;
     default:
       // Unknown type, error and exit
       logError("found unknown event type %d while reading event file\n",
@@ -159,6 +185,10 @@ const char *eventTypeToString(event_type_t type) {
       return "fert";
     case TILLAGE:
       return "till";
+    case LEAFON:
+      return "leafon";
+    case LEAFOFF:
+      return "leafoff";
     default:
       logError("unknown event type in eventTypeToString (%d)", type);
       exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM);
@@ -168,15 +198,26 @@ const char *eventTypeToString(event_type_t type) {
 event_type_t eventStringToType(const char *eventTypeStr) {
   if (strcmp(eventTypeStr, "irrig") == 0) {
     return IRRIGATION;
-  } else if (strcmp(eventTypeStr, "fert") == 0) {
+  }
+  if (strcmp(eventTypeStr, "fert") == 0) {
     return FERTILIZATION;
-  } else if (strcmp(eventTypeStr, "plant") == 0) {
+  }
+  if (strcmp(eventTypeStr, "plant") == 0) {
     return PLANTING;
-  } else if (strcmp(eventTypeStr, "till") == 0) {
+  }
+  if (strcmp(eventTypeStr, "till") == 0) {
     return TILLAGE;
-  } else if (strcmp(eventTypeStr, "harv") == 0) {
+  }
+  if (strcmp(eventTypeStr, "harv") == 0) {
     return HARVEST;
   }
+  if (strcmp(eventTypeStr, "leafon") == 0) {
+    return LEAFON;
+  }
+  if (strcmp(eventTypeStr, "leafoff") == 0) {
+    return LEAFOFF;
+  }
+
   return UNKNOWN_EVENT;
 }
 
@@ -193,6 +234,18 @@ static void checkEventLineTruncation(const char *line, size_t len) {
   }
 }
 
+void checkForCalculatedLeafEvents(void) {
+  // We have a leaf event in events.in, so make sure we are not also
+  // calculating leaf events
+  if (ctx.gdd || ctx.soilPhenol || params.leafOnDay > 0 ||
+      params.leafOffDay > 0) {
+    logError("calculated leaf events (via leafOnDay/leafOffDay params or "
+             "gdd/soil-phenol command-line options) are not compatible "
+             "with user-specified leaf events in event file\n");
+    exit(EXIT_CODE_BAD_PARAMETER_VALUE);
+  }
+}
+
 EventNode *readEventData(const char *eventFile) {
   int year, day, eventType;
   int currYear, currDay;
@@ -202,6 +255,7 @@ EventNode *readEventData(const char *eventFile) {
   char line[EVENT_LINE_SIZE];
   EventNode *curr, *next;
   EventNode *newEvents = NULL;
+  int hasCheckedLeafEvents = 0;
 
   // Check for a non-empty file
   if (access(eventFile, F_OK) != 0) {
@@ -238,6 +292,13 @@ EventNode *readEventData(const char *eventFile) {
     exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM);
   }
 
+  if (eventType == LEAFOFF || eventType == LEAFON) {
+    // If we have a leaf event, make sure we aren't also looking for
+    // calculated leaf events.
+    checkForCalculatedLeafEvents();
+    hasCheckedLeafEvents = 1;
+  }
+
   newEvents = createEventNode(year, day, eventType, eventParamsStr);
   next = newEvents;
   currYear = year;
@@ -264,12 +325,20 @@ EventNode *readEventData(const char *eventFile) {
       exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM);
     }
 
+    if (eventType == LEAFOFF || eventType == LEAFON) {
+      // If we have a leaf event, make sure we aren't also looking for
+      // calculated leaf events.
+      if (!hasCheckedLeafEvents) {
+        checkForCalculatedLeafEvents();
+        hasCheckedLeafEvents = 1;
+      }
+    }
+
     if ((year < currYear) || ((year == currYear) && (day < currDay))) {
-      // clang-format off
       logError("reading event file: last event was at (%d, %d), next event is "
-               "at (%d, %d)\n", currYear, currDay, year, day);
+               "at (%d, %d)\n",
+               currYear, currDay, year, day);
       logError("event records must be in time-ascending order\n");
-      // clang-format on
       exit(EXIT_CODE_INPUT_FILE_ERROR);
     }
 
@@ -363,30 +432,9 @@ int isFirstEventBefore(int year, int day) {
   return firstEvent->day < day;
 }
 
-void resetEventFluxes(void) {
-  fluxes.eventLeafC = 0.0;
-  fluxes.eventWoodC = 0.0;
-  fluxes.eventFineRootC = 0.0;
-  fluxes.eventCoarseRootC = 0.0;
-  fluxes.eventEvap = 0.0;
-  fluxes.eventSoilWater = 0.0;
-  fluxes.eventSoilC = 0.0;
-  fluxes.eventLitterC = 0.0;
-  fluxes.eventMinN = 0.0;
-  fluxes.eventSoilOrgN = 0.0;
-  fluxes.eventLitterN = 0.0;
-
-  // mass balance
-  fluxes.eventInputC = 0.0;
-  fluxes.eventOutputC = 0.0;
-  fluxes.eventInputN = 0.0;
-  fluxes.eventOutputN = 0.0;
-}
-
 void processEvents(void) {
-  // Set all event fluxes to zero, as these have no memory from one step to
-  // the next
-  resetEventFluxes();
+  // Event fluxes have all been reset to zero at the start of the time step,
+  // so we can just add to them as needed
 
   // If event starts off NULL, this function will just fall through, as it
   // should.
@@ -537,15 +585,20 @@ void processEvents(void) {
                          fracRB);
           fluxes.eventOutputN += outputN / climLen;
         }
+        // clang-format off
         writeEventOut(
-            gEvent, 10, "fluxes.eventSoilC", soilAdd / climLen,
-            "fluxes.eventLitterC", litterAdd / climLen, "fluxes.eventLeafC",
-            leafDelta / climLen, "fluxes.eventWoodC", woodDelta / climLen,
+            gEvent, 10,
+            "fluxes.eventSoilC", soilAdd / climLen,
+            "fluxes.eventLitterC", litterAdd / climLen,
+            "fluxes.eventLeafC", leafDelta / climLen,
+            "fluxes.eventWoodC", woodDelta / climLen,
             "fluxes.eventFineRootC", fineDelta / climLen,
             "fluxes.eventCoarseRootC", coarseDelta / climLen,
-            "fluxes.eventSoilOrgN", soilNAdd / climLen, "fluxes.eventLitterN",
-            litterNAdd / climLen, "fluxes.eventOutputC", outputC / climLen,
+            "fluxes.eventSoilOrgN", soilNAdd / climLen,
+            "fluxes.eventLitterN", litterNAdd / climLen,
+            "fluxes.eventOutputC", outputC / climLen,
             "fluxes.eventOutputN", outputN / climLen);
+        // clang-format on
       } break;
       case TILLAGE: {
         // BIG NOTE: this is the one event type that is NOT modeled as a flux;
@@ -561,8 +614,12 @@ void processEvents(void) {
       case FERTILIZATION: {
         const FertilizationParams *fertParams = gEvent->eventParams;
         const double orgC = fertParams->orgC;
-        double orgN = fertParams->orgN;
-        double minN = fertParams->minN;
+        double orgN = 0.0;
+        double minN = 0.0;
+        if (ctx.nitrogenCycle) {
+          orgN = fertParams->orgN;
+          minN = fertParams->minN;
+        }
         if (ctx.litterPool) {
           fluxes.eventLitterC += orgC / climLen;
         } else {
@@ -583,10 +640,47 @@ void processEvents(void) {
           fluxes.eventInputN += (orgN + minN) / climLen;
         }
 
-        writeEventOut(gEvent, 5, "fluxes.eventOrgN", orgN / climLen,
-                      "fluxes.eventLitterC", orgC / climLen, "fluxes.eventMinN",
-                      minN / climLen, "fluxes.eventInputC", orgC / climLen,
-                      "fluxes.eventInputN", (orgN + minN) / climLen);
+        // clang-format off
+        writeEventOut(gEvent, 6,
+          "fluxes.eventLitterC", ctx.litterPool ? orgC / climLen : 0.0,
+          "fluxes.eventSoilC", ctx.litterPool ? 0.0 : orgC / climLen,
+          "fluxes.eventMinN", minN / climLen,
+          "fluxes.eventLitterN", orgN / climLen,
+          "fluxes.eventInputC", orgC / climLen,
+          "fluxes.eventInputN", (orgN + minN) / climLen);
+        // clang-format on
+      } break;
+      case LEAFON: {
+        double leafOn = params.leafGrowth;
+        fluxes.eventLeafOnCreation += leafOn / climLen;
+
+        // Nitrogen is handled implicitly by relative CN ratios. Missing N
+        // from low-N wood to higher-N leaves is accounted for in
+        // calcNFixationAndUptakeFluxes() via calcPlantNDemand()
+
+        // Unlike planting, this is NOT a system input, so no adjustments to
+        // eventInputC or eventInputN
+
+        // ALSO, unlike all other events, we don't write the event here, as N
+        // limitation may change the amount
+      } break;
+      case LEAFOFF: {
+        double leafOff = envi.plantLeafC * params.fracLeafFall;
+        fluxes.eventLeafOffLitter += leafOff / climLen;
+
+        double litterNAdd = 0.0;
+        if (ctx.nitrogenCycle) {
+          // Nitrogen - need to account for leaf N moving to litter, as with
+          // harvests
+          litterNAdd = leafOff / params.leafCN;
+          fluxes.eventLitterN += litterNAdd / climLen;
+        }
+
+        // clang-format off
+        writeEventOut(gEvent, 2,
+          "fluxes.eventLeafOffLitter", leafOff / climLen,
+          "fluxes.eventLitterN", litterNAdd / climLen);
+        // clang-format on
       } break;
       default:
         logError("Unknown event type (%d) in processEvents()\n", gEvent->type);
@@ -609,6 +703,16 @@ void updatePoolsForEvents(void) {
     envi.litterC += fluxes.eventLitterC * climate->length;
   }
 
+  // Leaf on and off events
+  envi.plantWoodC -= fluxes.eventLeafOnCreation * climate->length;
+  envi.plantLeafC += (fluxes.eventLeafOnCreation - fluxes.eventLeafOffLitter) *
+                     climate->length;
+  if (ctx.litterPool) {
+    envi.litterC += fluxes.eventLeafOffLitter * climate->length;
+  } else {
+    envi.soilC += fluxes.eventLeafOffLitter * climate->length;
+  }
+
   // Harvest and planting events
   envi.coarseRootC += fluxes.eventCoarseRootC * climate->length;
   envi.fineRootC += fluxes.eventFineRootC * climate->length;