diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 5fe7e618..c7e20d6b 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -26,6 +26,8 @@ sections to include in release notes:
 ### Added
 
 ### Fixed
+- Terminating harvest events will now remove all biomass at end of time step, instead of allowing small amounts to remain in the plant pools due to growth in that step. (#400)
+- Leaf-off mechanics refined to ensure that leaf-off events do not remove more leaf carbon than is available in the leaf pool. (#400)
 
 ### Changed
 - SIPNET will now error instead of warning when an environment pool goes negative. This is a change from previous behavior where SIPNET would log a warning and continue running. (#398)
diff --git a/src/common/exitCodes.h b/src/common/exitCodes.h
index 1f17776b..2c399490 100644
--- a/src/common/exitCodes.h
+++ b/src/common/exitCodes.h
@@ -23,7 +23,8 @@ typedef enum {
   EXIT_CODE_FILE_OPEN_OR_READ_ERROR = 6,
   EXIT_CODE_INTERNAL_ERROR = 7,
   EXIT_CODE_BAD_CLI_ARGUMENT = 8,
-  EXIT_CODE_BAD_RESTART_PARAMETER = 9
+  EXIT_CODE_BAD_RESTART_PARAMETER = 9,
+  EXIT_CODE_MEMORY_ALLOCATION_FAILURE = 10
 } exit_code_t;
 
 #endif
diff --git a/src/common/util.c b/src/common/util.c
index fb087421..632c8e75 100644
--- a/src/common/util.c
+++ b/src/common/util.c
@@ -76,3 +76,109 @@ double calcRatio(const double num, const double den) {
 
 // For global linkage
 extern inline double unitClip(double preClip);
+
+DynamicString *dsCreate(size_t initial_capacity) {
+  DynamicString *ds = malloc(sizeof(DynamicString));
+  if (!ds)
+    return NULL;
+
+  // Ensure we have room for at least a null terminator
+  ds->capacity = (initial_capacity > 0) ? initial_capacity : 16;
+  ds->buffer = malloc(ds->capacity * sizeof(char));
+
+  if (!ds->buffer) {
+    free(ds);
+    return NULL;
+  }
+
+  ds->buffer[0] = '\0';  // Start with an empty string
+  ds->length = 0;
+  return ds;
+}
+
+int dsAppend(DynamicString *ds, const char *str) {
+  if (!ds || !str)
+    return 0;
+
+  size_t append_len = strlen(str);
+  // +1 is crucial to ensure room for the null terminator
+  size_t needed_capacity = ds->length + append_len + 1;
+
+  // Double capacity until it's large enough for the new content
+  if (needed_capacity > ds->capacity) {
+    size_t new_capacity = ds->capacity * 2;
+    while (new_capacity < needed_capacity) {
+      new_capacity *= 2;
+    }
+
+    // Safely reallocate using a temporary pointer
+    char *temp = realloc(ds->buffer, new_capacity * sizeof(char));
+    if (!temp) {
+      return 0;  // Reallocation failed, original data remains intact
+    }
+
+    ds->buffer = temp;
+    ds->capacity = new_capacity;
+  }
+
+  // Copy the new string over the old null terminator
+  memcpy(ds->buffer + ds->length, str, append_len);
+  ds->length += append_len;
+  ds->buffer[ds->length] = '\0';  // Manually place the new null terminator
+
+  return 1;  // Success
+}
+
+int dsAppendFormatted(DynamicString *ds, const char *format, ...) {
+  if (!ds || !format)
+    return 0;
+
+  // 1. Determine how much space the formatted text needs
+  va_list args;
+  va_start(args, format);
+  // Make a copy of args because vsnprintf consumes the list
+  va_list args_copy;
+  va_copy(args_copy, args);
+
+  // Pass NULL and 0 to just count the characters needed
+  int formatted_len = vsnprintf(NULL, 0, format, args_copy);
+  va_end(args_copy);
+
+  if (formatted_len < 0) {
+    va_end(args);
+    return 0;  // Formatting error occurred
+  }
+
+  // 2. Ensure the buffer is big enough
+  size_t needed_capacity = ds->length + (size_t)formatted_len + 1;
+  if (needed_capacity > ds->capacity) {
+    size_t new_capacity = ds->capacity * 2;
+    while (new_capacity < needed_capacity) {
+      new_capacity *= 2;
+    }
+
+    char *temp = realloc(ds->buffer, new_capacity * sizeof(char));
+    if (!temp) {
+      va_end(args);
+      return 0;  // Reallocation failed
+    }
+    ds->buffer = temp;
+    ds->capacity = new_capacity;
+  }
+
+  // 3. Write the formatted string directly into the builder's buffer
+  // Write directly to the position of the current null terminator
+  vsnprintf(ds->buffer + ds->length, (size_t)formatted_len + 1, format, args);
+  va_end(args);
+
+  // 4. Update the length of the string builder
+  ds->length += (size_t)formatted_len;
+  return 1;
+}
+
+void dsFree(DynamicString *ds) {
+  if (ds) {
+    free(ds->buffer);
+    free(ds);
+  }
+}
diff --git a/src/common/util.h b/src/common/util.h
index c8a31f26..8344d7a3 100644
--- a/src/common/util.h
+++ b/src/common/util.h
@@ -9,7 +9,10 @@
 #define UTIL_H
 
 #include <math.h>
+#include <stdarg.h>  // Required for va_list, va_start, va_end
 #include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 
 #define TINY 0.000001  // to avoid those nasty divide-by-zero errors
 
@@ -37,4 +40,23 @@ double calcRatio(double num, double den);
  */
 inline double unitClip(double preClip) { return fmin(fmax(preClip, 0.0), 1.0); }
 
+// DYNAMIC STRING
+typedef struct DynamicStringStruct {
+  char *buffer;  // Pointer to the character array
+  size_t length;  // Number of characters currently in the string
+  size_t capacity;  // Total allocated space (including null terminator)
+} DynamicString;
+
+// Initialize the builder with a reasonable default capacity
+DynamicString *dsCreate(size_t initial_capacity);
+
+// Append text, automatically doubling capacity if needed
+int dsAppend(DynamicString *ds, const char *str);
+
+// Append formatted text using printf-style syntax
+int dsAppendFormatted(DynamicString *ds, const char *format, ...);
+
+// Free all memory associated with the builder
+void dsFree(DynamicString *ds);
+
 #endif
diff --git a/src/sipnet/balance.c b/src/sipnet/balance.c
index 1c43608c..2efe55ce 100644
--- a/src/sipnet/balance.c
+++ b/src/sipnet/balance.c
@@ -47,7 +47,8 @@ void updateBalanceTrackerPostClamp(void) {
   balanceTracker.clampedC = balanceTracker.finalC - balanceTracker.postTotalC;
   if (balanceTracker.clampedC < -EPS) {
     // This shouldn't happen, by construction
-    logInternalError("Non-negative clamping has cause carbon loss\n");
+    logInternalError("Non-negative clamping has cause carbon loss %f\n",
+                     balanceTracker.clampedC);
   }
   if (balanceTracker.clampedC < EPS) {
     balanceTracker.clampedC = 0;
@@ -56,7 +57,8 @@ void updateBalanceTrackerPostClamp(void) {
   balanceTracker.clampedN = balanceTracker.finalN - balanceTracker.postTotalN;
   if (balanceTracker.clampedN < -EPS) {
     // This shouldn't happen, by construction
-    logInternalError("Non-negative clamping has cause nitrogen loss\n");
+    logInternalError("Non-negative clamping has cause nitrogen loss %f\n",
+                     balanceTracker.clampedN);
   }
   if (balanceTracker.clampedN < EPS) {
     balanceTracker.clampedN = 0;
diff --git a/src/sipnet/debug_log.c b/src/sipnet/debug_log.c
index fa8e99ba..6f3119a8 100644
--- a/src/sipnet/debug_log.c
+++ b/src/sipnet/debug_log.c
@@ -4,11 +4,44 @@
 
 #include "debug_log.h"
 
+#include "events.h"
 #include "common/exitCodes.h"
 #include "common/logging.h"
 #include "common/context.h"
 #include "common/util.h"
 #include "state.h"
+
+#define NUM_LOGGED_ENVI_FIELDS 13
+#define NUM_LOGGED_FLUX_FIELDS 59
+#define NUM_LOGGED_TRACKER_FIELDS 33
+#define NUM_LOGGED_PHEN_TRACKER_FIELDS 3
+#define NUM_LOGGED_SURVIVAL_FIELDS 1
+#define NUM_LOGGED_EVENT_TRACKER_FIELDS 7
+
+#define DEBUG_LAYOUT_ENVI_SIZE (8 * NUM_LOGGED_ENVI_FIELDS)
+#define DEBUG_LAYOUT_FLUX_SIZE (8 * NUM_LOGGED_FLUX_FIELDS)
+// The Trackers struct is not all doubles, but the int(s) get padded to 8 bytes
+#define DEBUG_LAYOUT_TRACKER_SIZE (8 * NUM_LOGGED_TRACKER_FIELDS)
+#define DEBUG_LAYOUT_PHEN_SIZE (4 * NUM_LOGGED_PHEN_TRACKER_FIELDS)
+#define DEBUG_LAYOUT_SURVIVAL_SIZE (4 * NUM_LOGGED_SURVIVAL_FIELDS)
+#define DEBUG_LAYOUT_EVENT_SIZE (8 * NUM_LOGGED_EVENT_TRACKER_FIELDS)
+
+_Static_assert(sizeof(Envi) == DEBUG_LAYOUT_ENVI_SIZE,
+               "Debug log schema drift: Envi changed; update debug_log.c");
+_Static_assert(sizeof(Fluxes) == DEBUG_LAYOUT_FLUX_SIZE,
+               "Debug log schema drift: Fluxes changed; update debug_log.c");
+_Static_assert(sizeof(Trackers) == DEBUG_LAYOUT_TRACKER_SIZE,
+               "Debug log schema drift: Trackers changed; update debug_log.c");
+_Static_assert(
+    sizeof(PhenologyTrackers) == DEBUG_LAYOUT_PHEN_SIZE,
+    "Debug log schema drift: PhenologyTrackers changed; update debug_log.c");
+_Static_assert(
+    sizeof(PlantSurvivalTracker) == DEBUG_LAYOUT_SURVIVAL_SIZE,
+    "Debug log schema drift: PlantSurvival changed; update debug_log.c");
+_Static_assert(
+    sizeof(EventTrackers) == DEBUG_LAYOUT_EVENT_SIZE,
+    "Debug log schema drift: EventTrackers changed; update debug_log.c");
+
 typedef enum DebugFieldType {
   DEBUG_FIELD_INT = 0,
   DEBUG_FIELD_DOUBLE = 1
@@ -20,18 +53,13 @@ typedef struct DebugField {
   const void *value;
 } DebugField;
 
-#define NUM_LOGGED_ENVI_FIELDS 13
-#define NUM_LOGGED_FLUX_FIELDS 56
-#define NUM_LOGGED_TRACKER_FIELDS 33
-#define NUM_LOGGED_PHEN_TRACKER_FIELDS 3
-#define NUM_LOGGED_SURVIVAL_FIELDS 1
-
 typedef struct DebugFieldArrays {
   DebugField enviDF[NUM_LOGGED_ENVI_FIELDS];
   DebugField fluxDF[NUM_LOGGED_FLUX_FIELDS];
   DebugField trackerDF[NUM_LOGGED_TRACKER_FIELDS];
   DebugField phenoDF[NUM_LOGGED_PHEN_TRACKER_FIELDS];
   DebugField survivalDF[NUM_LOGGED_SURVIVAL_FIELDS];
+  DebugField eventDF[NUM_LOGGED_EVENT_TRACKER_FIELDS];
 } DebugFieldArrays;
 
 static DebugFieldArrays *debugFields = NULL;
@@ -61,7 +89,11 @@ void initDebugArrays() {
   debugFields->enviDF[ind++] = (DebugField){"soilOrgN", DEBUG_FIELD_DOUBLE, &envi.soilOrgN},
   debugFields->enviDF[ind++] = (DebugField){"litterN", DEBUG_FIELD_DOUBLE, &envi.litterN},
   debugFields->enviDF[ind++] = (DebugField){"plantStorageN", DEBUG_FIELD_DOUBLE, &envi.plantStorageN},
-  debugFields->enviDF[ind  ] = (DebugField){"plantCAccountingDelta", DEBUG_FIELD_DOUBLE,&envi.plantCAccountingDelta};
+  debugFields->enviDF[ind++] = (DebugField){"plantCAccountingDelta", DEBUG_FIELD_DOUBLE,&envi.plantCAccountingDelta};
+  if (ind != NUM_LOGGED_ENVI_FIELDS) {
+    logInternalError("Debug log array size mismatch: enviDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
   debugFields->fluxDF[ind++] = (DebugField){"photosynthesis", DEBUG_FIELD_DOUBLE, &fluxes.photosynthesis};
@@ -90,6 +122,7 @@ void initDebugArrays() {
   debugFields->fluxDF[ind++] = (DebugField){"woodCreation", DEBUG_FIELD_DOUBLE, &fluxes.woodCreation},
   debugFields->fluxDF[ind++] = (DebugField){"leafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.leafOnCreation},
   debugFields->fluxDF[ind++] = (DebugField){"leafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.leafOnCreationFromWood},
+  debugFields->fluxDF[ind++] = (DebugField){"leafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.leafOffLitter},
   debugFields->fluxDF[ind++] = (DebugField){"nVolatilization", DEBUG_FIELD_DOUBLE, &fluxes.nVolatilization},
   debugFields->fluxDF[ind++] = (DebugField){"nLeaching", DEBUG_FIELD_DOUBLE, &fluxes.nLeaching},
   debugFields->fluxDF[ind++] = (DebugField){"nOrgSoil", DEBUG_FIELD_DOUBLE, &fluxes.nOrgSoil},
@@ -101,6 +134,7 @@ void initDebugArrays() {
   debugFields->fluxDF[ind++] = (DebugField){"reductionNResorption", DEBUG_FIELD_DOUBLE, &fluxes.reductionNResorption},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafC", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafC},
   debugFields->fluxDF[ind++] = (DebugField){"eventWoodC", DEBUG_FIELD_DOUBLE, &fluxes.eventWoodC},
+  debugFields->fluxDF[ind++] = (DebugField){"eventAccountingC", DEBUG_FIELD_DOUBLE, &fluxes.eventAccountingC},
   debugFields->fluxDF[ind++] = (DebugField){"eventFineRootC", DEBUG_FIELD_DOUBLE, &fluxes.eventFineRootC},
   debugFields->fluxDF[ind++] = (DebugField){"eventCoarseRootC", DEBUG_FIELD_DOUBLE, &fluxes.eventCoarseRootC},
   debugFields->fluxDF[ind++] = (DebugField){"eventEvap", DEBUG_FIELD_DOUBLE, &fluxes.eventEvap},
@@ -116,10 +150,15 @@ void initDebugArrays() {
   debugFields->fluxDF[ind++] = (DebugField){"eventOutputN", DEBUG_FIELD_DOUBLE, &fluxes.eventOutputN},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreation", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreation},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafOnCreationFromWood", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOnCreationFromWood},
-  debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitter", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitter},
+  debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitterC", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterC},
+  debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffLitterN", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffLitterN},
   debugFields->fluxDF[ind++] = (DebugField){"eventLeafOffNResorption", DEBUG_FIELD_DOUBLE, &fluxes.eventLeafOffNResorption},
   debugFields->fluxDF[ind++] = (DebugField){"soilMethane", DEBUG_FIELD_DOUBLE, &fluxes.soilMethane},
-  debugFields->fluxDF[ind  ] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane};
+  debugFields->fluxDF[ind++] = (DebugField){"litterMethane", DEBUG_FIELD_DOUBLE, &fluxes.litterMethane};
+  if (ind != NUM_LOGGED_FLUX_FIELDS) {
+    logInternalError("Debug log array size mismatch: fluxDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
   debugFields->trackerDF[ind++] = (DebugField){"gpp", DEBUG_FIELD_DOUBLE, &trackers.gpp},
@@ -154,15 +193,40 @@ void initDebugArrays() {
   debugFields->trackerDF[ind++] = (DebugField){"nLeaching", DEBUG_FIELD_DOUBLE, &trackers.nLeaching},
   debugFields->trackerDF[ind++] = (DebugField){"nFixation", DEBUG_FIELD_DOUBLE, &trackers.nFixation},
   debugFields->trackerDF[ind++] = (DebugField){"nUptake", DEBUG_FIELD_DOUBLE, &trackers.nUptake};
-  debugFields->trackerDF[ind  ] = (DebugField){"meanNPP", DEBUG_FIELD_DOUBLE, &trackers.meanNPP};
+  debugFields->trackerDF[ind++] = (DebugField){"meanNPP", DEBUG_FIELD_DOUBLE, &trackers.meanNPP};
+  if (ind != NUM_LOGGED_TRACKER_FIELDS) {
+    logInternalError("Debug log array size mismatch: trackerDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
   debugFields->phenoDF[ind++] = (DebugField){"didLeafGrowth", DEBUG_FIELD_INT, &phenologyTrackers.didLeafGrowth},
   debugFields->phenoDF[ind++] = (DebugField){"didLeafFall", DEBUG_FIELD_INT, &phenologyTrackers.didLeafFall},
-  debugFields->phenoDF[ind  ] = (DebugField){"lastYear", DEBUG_FIELD_INT, &phenologyTrackers.lastYear};
+  debugFields->phenoDF[ind++] = (DebugField){"lastYear", DEBUG_FIELD_INT, &phenologyTrackers.lastYear};
+  if (ind != NUM_LOGGED_PHEN_TRACKER_FIELDS) {
+    logInternalError("Debug log array size mismatch: phenoDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
 
   ind = 0;
-  debugFields->survivalDF[ind] = (DebugField){"isAlive", DEBUG_FIELD_INT, &plantSurvivalTracker.isAlive};
+  debugFields->survivalDF[ind++] = (DebugField){"isAlive", DEBUG_FIELD_INT, &plantSurvivalTracker.isAlive};
+  if (ind != NUM_LOGGED_SURVIVAL_FIELDS) {
+    logInternalError("Debug log array size mismatch: survivalDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
+
+  ind = 0;
+  debugFields->eventDF[ind++] = (DebugField){"d_till_mod", DEBUG_FIELD_DOUBLE, &eventTrackers.d_till_mod};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemoved", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemoved};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferred", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferred};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemovedAbove", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedAbove};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracRemovedBelow", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracRemovedBelow};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferredAbove", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredAbove};
+  debugFields->eventDF[ind++] = (DebugField){"ht.totalFracTransferredBelow", DEBUG_FIELD_DOUBLE, &eventTrackers.harvestTrackers.totalFracTransferredBelow};
+  if (ind != NUM_LOGGED_EVENT_TRACKER_FIELDS) {
+    logInternalError("Debug log array size mismatch: eventDF\n");
+    exit(EXIT_CODE_INTERNAL_ERROR);
+  }
   // clang-format on
 }
 
@@ -278,6 +342,8 @@ void outputDebugHeaders(DebugLogFiles *debugLogFiles) {
     outputDebugFieldHeader(debugLogFiles->trackers, "s.",
                            debugFields->survivalDF, NUM_LOGGED_SURVIVAL_FIELDS,
                            0);
+    outputDebugFieldHeader(debugLogFiles->trackers, "et.", debugFields->eventDF,
+                           NUM_LOGGED_EVENT_TRACKER_FIELDS, 0);
     fprintf(debugLogFiles->trackers, "\n");
   }
 }
@@ -307,6 +373,9 @@ void outputDebugState(DebugLogFiles *debugLogFiles, int year, int day,
     outputDebugFieldValues(debugLogFiles->trackers, year, day, time,
                            debugFields->survivalDF, NUM_LOGGED_SURVIVAL_FIELDS,
                            0);
+    outputDebugFieldValues(debugLogFiles->trackers, year, day, time,
+                           debugFields->eventDF,
+                           NUM_LOGGED_EVENT_TRACKER_FIELDS, 0);
     fprintf(debugLogFiles->trackers, "\n");
   }
 }
diff --git a/src/sipnet/events.c b/src/sipnet/events.c
index cc99ce0b..6b7f9765 100644
--- a/src/sipnet/events.c
+++ b/src/sipnet/events.c
@@ -43,6 +43,8 @@ EventNode *createEventNode(int year, int day, int eventType,
   newEvent->year = year;
   newEvent->day = day;
   newEvent->type = eventType;
+  newEvent->numLogParamPairs = 0;
+  newEvent->logLine = dsCreate(0);
 
   switch (eventType) {
     case HARVEST: {
@@ -58,7 +60,13 @@ EventNode *createEventNode(int year, int day, int eventType,
       // Validate the params
       if ((fracRA + fracTA > 1) || (fracRB + fracTB > 1)) {
         logError("invalid harvest newEvent for year %d day %d; above and below "
-                 "must each add to 1 or less",
+                 "must each add to 1 or less\n",
+                 year, day);
+        exit(EXIT_CODE_BAD_PARAMETER_VALUE);
+      }
+      if (fracRA < 0.0 || fracRB < 0.0 || fracTA < 0.0 || fracTB < 0.0) {
+        logError("invalid harvest newEvent for year %d day %d; fractions must "
+                 "be non-negative\n",
                  year, day);
         exit(EXIT_CODE_BAD_PARAMETER_VALUE);
       }
@@ -369,51 +377,64 @@ EventNode *readEventData(const char *eventFile) {
 void openEventOutFile(const char *eventOutFilePath, int printHeader) {
   eventOutFile = openFile(eventOutFilePath, "w");
   if (printHeader) {
-    // Use format string analogous to the one in writeEventOut for
+    // Use format string analogous to the one in writeEventsOut for
     // better alignment (won't be perfect, but definitely better)
     fprintf(eventOutFile, "%4s  %3s  %-7s  %s", "year", "day", "type",
             "param_name=delta[,param_name=delta,...]\n");
   }
 }
 
-void doWriteEventOut(int year, int day, const char *type, int numParams,
-                     va_list args) {
-  int ind = 0;
-
-  // Spec:
-  // year day event_type <param name=delta>[,<param name>=<delta>,...]
-
-  // Standard prefix for all
-  fprintf(eventOutFile, "%4d  %3d  %-7s  ", year, day, type);
-
-  // For debugging on linux
-  char *param;
-  double val;
-  // Variable output per oneEvent type
-  for (ind = 0; ind < numParams - 1; ind++) {
-    // For debugging on linux
-    param = va_arg(args, char *);
-    val = va_arg(args, double);
-    fprintf(eventOutFile, "%s=%-.2f,", param, val);
-  }
-  param = va_arg(args, char *);
-  val = va_arg(args, double);
-  fprintf(eventOutFile, "%s=%-.2f\n", param, val);
-}
-
-void writeEventOut(EventNode *oneEvent, int numParams, ...) {
+void appendLog(EventNode *event, int numParams, ...) {
   va_list args;
   va_start(args, numParams);
-  doWriteEventOut(oneEvent->year, oneEvent->day,
-                  eventTypeToString(oneEvent->type), numParams, args);
+  for (int ind = 0; ind < numParams; ind++) {
+    char *param = va_arg(args, char *);
+    double val = va_arg(args, double);
+    int success = dsAppendFormatted(event->logLine, "%s=%-.2f,", param, val);
+    if (!success) {
+      logError("appending event log line failed\n");
+      exit(EXIT_CODE_INTERNAL_ERROR);
+    }
+  }
   va_end(args);
 }
 
+void writeEventsOut(void) {
+  // We move the global events pointer here (instead of a copy as in the
+  // processEvents functions), as this is the last pass
+  const int climYear = climate->year;
+  const int climDay = climate->day;
+  while (gEvent != NULL && gEvent->year <= climYear && gEvent->day <= climDay) {
+    // Change last char to a newline if it is a comma
+    // Get the current length of the string
+    char *log = gEvent->logLine->buffer;
+    size_t len = strlen(log);
+    if (len > 0 && log[len - 1] == ',') {
+      log[len - 1] = '\n';
+    } else {
+      dsAppend(gEvent->logLine, "\n");
+    }
+    fprintf(eventOutFile, "%4d  %3d  %-7s  %s", gEvent->year, gEvent->day,
+            eventTypeToString(gEvent->type), log);
+    gEvent = gEvent->nextEvent;
+  }
+}
+
 void writeComputedEventOut(int year, int day, const char *type, int numParams,
                            ...) {
   va_list args;
   va_start(args, numParams);
-  doWriteEventOut(year, day, type, numParams, args);
+  // Standard prefix for all
+  fprintf(eventOutFile, "%4d  %3d  %-7s  ", year, day, type);
+
+  // Variable output per oneEvent type
+  for (int ind = 0; ind < numParams; ind++) {
+    char *param = va_arg(args, char *);
+    double val = va_arg(args, double);
+    char suffix = (ind == numParams - 1) ? '\n' : ',';
+    fprintf(eventOutFile, "%s=%-.2f%c", param, val, suffix);
+  }
+
   va_end(args);
 }
 
@@ -446,7 +467,9 @@ int isFirstEventBefore(int year, int day) {
   return firstEvent->day < day;
 }
 
-void processEvents(void) {
+EventNode *getCurrentEvent(void) { return gEvent; }
+
+void processEventsForCarbon(EventNode *event) {
   // Event fluxes have all been reset to zero at the start of the time step,
   // so we can just add to them as needed
 
@@ -465,24 +488,27 @@ void processEvents(void) {
   }
 
   // Reset harvest tracking
-  eventTrackers.harvestFracRemoved = 0;
-  eventTrackers.harvestFracTransferred = 0;
+  eventTrackers.harvestTrackers = (HarvestTrackers){0};
 
-  while (gEvent != NULL && gEvent->year <= climYear && gEvent->day <= climDay) {
+  // Make sure we don't have more than 100% leaves fall for leaf-off
+  double totalLeafFallFrac = 0.0;
+
+  while (event != NULL && event->year <= climYear && event->day <= climDay) {
     // The events file has been tested on read, so we know this event list
     // should be in chrono order. However, we need to check to make sure the
     // current event is not in the past, as that would indicate an event that
     // did not have a corresponding climate file record.
-    if (gEvent->year < climYear || gEvent->day < climDay) {
+    if (event->year < climYear || event->day < climDay) {
       logError("Agronomic event found for year: %d day: %d that does not "
-               "have a corresponding record in the climate file\n",
-               gEvent->year, gEvent->day);
+               "have a corresponding record in the climate file "
+               "(current climate date is year %d day %d)\n",
+               event->year, event->day, climYear, climDay);
       exit(EXIT_CODE_INPUT_FILE_ERROR);
     }
 
-    switch (gEvent->type) {
+    switch (event->type) {
       case IRRIGATION: {
-        const IrrigationParams *irrParams = gEvent->eventParams;
+        const IrrigationParams *irrParams = event->eventParams;
         const double amount = irrParams->amountAdded;
         double soilAmount, evapAmount;
         if (irrParams->method == CANOPY) {
@@ -501,11 +527,13 @@ void processEvents(void) {
         }
         fluxes.eventEvap += evapAmount / climLen;
         fluxes.eventSoilWater += soilAmount / climLen;
-        writeEventOut(gEvent, 2, "eventSoilWater", soilAmount, "eventEvap",
-                      evapAmount);
+
+        appendLog(event, 2, "eventSoilWater", soilAmount, "eventEvap",
+                  evapAmount);
+
       } break;
       case PLANTING: {
-        const PlantingParams *plantParams = gEvent->eventParams;
+        const PlantingParams *plantParams = event->eventParams;
         const double leafC = plantParams->leafC;
         const double woodC = plantParams->woodC;
         const double fineRootC = plantParams->fineRootC;
@@ -522,31 +550,25 @@ void processEvents(void) {
 
         // MASS BALANCE: this is a system input
         const double inputC = leafC + woodC + fineRootC + coarseRootC;
-        double inputN = 0.0;
         fluxes.eventInputC += inputC / climLen;
-        if (ctx.nitrogenCycle) {
-          inputN = leafC / params.leafCN + woodC / params.woodCN +
-                   fineRootC / params.fineRootCN + coarseRootC / params.woodCN;
-          fluxes.eventInputN += inputN / climLen;
-        }
 
         // clang-format off
-        writeEventOut(gEvent, 6,
+        appendLog(event, 5,
                       "eventLeafC", leafC,
                       "eventWoodC", woodC,
                       "eventFineRootC", fineRootC,
                       "eventCoarseRootC", coarseRootC,
-                      "eventInputC", inputC,
-                      "eventInputN", inputN);
+                      "eventInputC", inputC);
         // clang-format on
+
       } break;
       case HARVEST: {
         // Harvest can both remove biomass and move biomass to the soil/litter
         // pools
-        const HarvestParams *harvParams = gEvent->eventParams;
+        const HarvestParams *harvParams = event->eventParams;
         const double fracRA = harvParams->fractionRemovedAbove;
-        const double fracTA = harvParams->fractionTransferredAbove;
         const double fracRB = harvParams->fractionRemovedBelow;
+        const double fracTA = harvParams->fractionTransferredAbove;
         const double fracTB = harvParams->fractionTransferredBelow;
         const double woodC = envi.plantWoodC + envi.plantCAccountingDelta;
 
@@ -554,11 +576,32 @@ void processEvents(void) {
         double aboveMass = woodC + envi.plantLeafC;
         double belowMass = envi.fineRootC + envi.coarseRootC;
         double totalMass = aboveMass + belowMass;
+        HarvestTrackers *ht = &eventTrackers.harvestTrackers;
         if (totalMass > TINY) {
           double massRemoved = fracRA * aboveMass + fracRB * belowMass;
           double massTransferred = fracTA * aboveMass + fracTB * belowMass;
-          eventTrackers.harvestFracRemoved += massRemoved / totalMass;
-          eventTrackers.harvestFracTransferred += massTransferred / totalMass;
+          ht->totalFracRemoved += massRemoved / totalMass;
+          ht->totalFracTransferred += massTransferred / totalMass;
+          ht->totalFracRemovedAbove += fracRA;
+          ht->totalFracRemovedBelow += fracRB;
+          ht->totalFracTransferredAbove += fracTA;
+          ht->totalFracTransferredBelow += fracTB;
+          if (ht->totalFracRemovedAbove + ht->totalFracTransferredAbove >
+                  1.0 + TINY ||
+              ht->totalFracRemovedBelow + ht->totalFracTransferredBelow >
+                  1.0 + TINY) {
+            logError("Harvest event(s) at year %d day %d has total above-ground"
+                     " or below-ground removal + transfer fraction > 1.0"
+                     " (above %.3f, below %.3f)\n",
+                     event->year, event->day,
+                     ht->totalFracRemovedAbove + ht->totalFracTransferredAbove,
+                     ht->totalFracRemovedBelow + ht->totalFracTransferredBelow);
+            exit(EXIT_CODE_BAD_PARAMETER_VALUE);
+          }
+        } else {
+          logWarning("Harvest event at year %d day %d has no biomass to remove "
+                     "or transfer\n",
+                     event->year, event->day);
         }
 
         // Litter increase
@@ -568,7 +611,9 @@ void processEvents(void) {
         // Pool reductions, counting both mass moved to litter and removed by
         // the harvest itself. Above-ground changes:
         const double leafDelta = -envi.plantLeafC * (fracRA + fracTA);
-        const double woodDelta = -woodC * (fracRA + fracTA);
+        const double woodDelta = -envi.plantWoodC * (fracRA + fracTA);
+        const double accountingDelta =
+            -envi.plantCAccountingDelta * (fracRA + fracTA);
         // Below-ground changes:
         const double fineDelta = -envi.fineRootC * (fracRB + fracTB);
         const double coarseDelta = -envi.coarseRootC * (fracRB + fracTB);
@@ -583,161 +628,266 @@ void processEvents(void) {
         fluxes.eventSoilC += soilAdd / climLen;
         fluxes.eventLeafC += leafDelta / climLen;
         fluxes.eventWoodC += woodDelta / climLen;
+        fluxes.eventAccountingC += accountingDelta / climLen;
         fluxes.eventFineRootC += fineDelta / climLen;
         fluxes.eventCoarseRootC += coarseDelta / climLen;
 
-        // No need to allocate to biomass N pools, we don't track that N
-        // explicitly. We do need to handle soil and litter N, though.
-        // Note: ctx.nitrogenCycle implies ctx.litterPool
-        // Litter N increase
-        double litterNAdd = 0.0;
-        double soilNAdd = 0.0;
-        if (ctx.nitrogenCycle) {
-          const double totalAbove = (envi.plantLeafC / params.leafCN) +
-                                    (envi.plantWoodC / params.woodCN);
-          const double totalBelow = (envi.fineRootC / params.fineRootCN) +
-                                    (envi.coarseRootC / params.woodCN);
-          litterNAdd = fracTA * totalAbove;
-          soilNAdd = fracTB * totalBelow;
-          fluxes.eventSoilOrgN += soilNAdd / climLen;
-          fluxes.eventLitterN += litterNAdd / climLen;
-        }
-
         // MASS BALANCE: removed fractions are system outputs
         const double outputC = ((woodC + envi.plantLeafC) * fracRA +
                                 (envi.fineRootC + envi.coarseRootC) * fracRB);
-        double outputN = 0.0;
         fluxes.eventOutputC += outputC / climLen;
-        if (ctx.nitrogenCycle) {
-          // just plantWoodC here, not woodC
-          outputN = (envi.plantWoodC / params.woodCN +
-                     envi.plantLeafC / params.leafCN) *
-                        fracRA +
-                    (envi.fineRootC / params.fineRootCN +
-                     envi.coarseRootC / params.woodCN) *
-                        fracRB;
-          fluxes.eventOutputN += outputN / climLen;
-        }
+
         // clang-format off
-        writeEventOut(
-            gEvent, 10,
+        appendLog(
+            event, 8,
             "eventSoilC", soilAdd,
             "eventLitterC", litterAdd,
             "eventLeafC", leafDelta,
             "eventWoodC", woodDelta,
+            "eventAccountingC", accountingDelta,
             "eventFineRootC", fineDelta,
             "eventCoarseRootC", coarseDelta,
-            "eventSoilOrgN", soilNAdd,
-            "eventLitterN", litterNAdd,
-            "eventOutputC", outputC,
-            "eventOutputN", outputN);
+            "eventOutputC", outputC);
         // clang-format on
+
       } break;
       case TILLAGE: {
         // BIG NOTE: this is the one event type that is NOT modeled as a flux;
         // see updateEventTrackers() for more
-        const TillageParams *tillParams = gEvent->eventParams;
+        const TillageParams *tillParams = event->eventParams;
         // Update the tillage mod for R_H calculations; this will be slowly
         // reduced by an exponential decay function. Note we add here, not set,
         // as there may be lingering effects from a prior tillage.
         eventTrackers.d_till_mod += tillParams->tillageEffect;
-        writeEventOut(gEvent, 1, "eventTrackers.d_till_mod",
-                      tillParams->tillageEffect);
+
+        appendLog(event, 1, "eventTrackers.d_till_mod",
+                  tillParams->tillageEffect);
+
       } break;
       case FERTILIZATION: {
-        const FertilizationParams *fertParams = gEvent->eventParams;
+        const FertilizationParams *fertParams = event->eventParams;
         const double orgC = fertParams->orgC;
-        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 {
           fluxes.eventSoilC += orgC / climLen;
         }
 
-        if (ctx.nitrogenCycle) {
-          // As the warning says in readEventData(), we ignore N when the
-          // nitrogen cycle model is off
-          // Implies ctx.litterPool
-          fluxes.eventLitterN += orgN / climLen;
-          fluxes.eventMinN += minN / climLen;
-        }
-
         // MASS BALANCE: this is a system input
         fluxes.eventInputC += orgC / climLen;
-        if (ctx.nitrogenCycle) {
-          fluxes.eventInputN += (orgN + minN) / climLen;
-        }
 
         // clang-format off
-        writeEventOut(gEvent, 6,
+        appendLog(event, 3,
           "eventLitterC", ctx.litterPool ? orgC : 0.0,
           "eventSoilC", ctx.litterPool ? 0.0 : orgC,
-          "eventMinN", minN,
-          "eventLitterN", orgN,
-          "eventInputC", orgC,
-          "eventInputN", (orgN + minN));
+          "eventInputC", orgC);
         // clang-format on
+
       } break;
       case LEAFON: {
         double leafOnFlux = params.leafGrowth / climLen;
+        double leafOnFluxFromWood = 0.0;
         checkLeafOnLimitation(&leafOnFlux);
         fluxes.eventLeafOnCreation += leafOnFlux;
         double totalSourceC = envi.plantWoodC + envi.coarseRootC;
         if (totalSourceC > TINY) {
-          fluxes.eventLeafOnCreationFromWood +=
-              leafOnFlux * envi.plantWoodC / totalSourceC;
+          leafOnFluxFromWood = leafOnFlux * envi.plantWoodC / totalSourceC;
+          fluxes.eventLeafOnCreationFromWood += leafOnFluxFromWood;
+        }
+        // Unlike planting, this is NOT a system input, so no adjustments to
+        // eventInputC
+
+        // clang-format off
+        appendLog(event, 2,
+          "eventLeafOnCreation", leafOnFlux * climLen,
+          "eventLeafOnCreationFromWood", leafOnFluxFromWood * climLen);
+        // clang-format on
+      } break;
+      case LEAFOFF: {
+        double frac = params.fracLeafFall;
+        totalLeafFallFrac += frac;
+        if (totalLeafFallFrac > 1.0) {
+          logError("Total leaf fall for leaf-off event(s) is greater than 100 "
+                   "percent (%.3f) for year %d day %d\n",
+                   totalLeafFallFrac, event->year, event->day);
+          exit(EXIT_CODE_BAD_PARAMETER_VALUE);
         }
 
+        double leafOff = envi.plantLeafC * params.fracLeafFall;
+        fluxes.eventLeafOffLitterC += leafOff / climLen;
+
+        appendLog(event, 1, "eventLeafOffLitter", leafOff);
+      } break;
+      case PLANTDEATH:
+        // There should be no way to get here, but covering our bases...
+        logWarning("PLANTDEATH event found for year %d day %d, but not "
+                   "implemented as an input event; ignoring\n",
+                   event->year, event->day);
+        break;
+      default:
+        logError("Unknown event type (%d) in processEvents()\n", event->type);
+        exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM);
+    }
+
+    event = event->nextEvent;
+  }
+}
+
+void processEventsForNitrogen(EventNode *event) {
+  if (!ctx.nitrogenCycle) {
+    return;
+  }
+
+  // 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.
+  const int climYear = climate->year;
+  const int climDay = climate->day;
+  const double climLen = climate->length;
+
+  // Checks performed in processEventsForCarbon are not repeated here unless
+  // necessary for nitrogen handling.
+
+  // Leaf-off events are a little special, as we only process the first one
+  // (see LEAFOFF below)
+  int firstLeafOff = 1;
+
+  while (event != NULL && event->year <= climYear && event->day <= climDay) {
+    switch (event->type) {
+      case IRRIGATION: {
+        // Nothing to do for irrigation events
+      } break;
+      case PLANTING: {
+        const PlantingParams *plantParams = event->eventParams;
+        const double leafC = plantParams->leafC;
+        const double woodC = plantParams->woodC;
+        const double fineRootC = plantParams->fineRootC;
+        const double coarseRootC = plantParams->coarseRootC;
+
+        // No need to allocate to biomass N pools, we don't track that N
+        // explicitly
+
+        // MASS BALANCE: this is a system input
+        double inputN = leafC / params.leafCN + woodC / params.woodCN +
+                        fineRootC / params.fineRootCN +
+                        coarseRootC / params.woodCN;
+        fluxes.eventInputN += inputN / climLen;
+
+        appendLog(event, 1, "eventInputN", inputN);
+      } break;
+      case HARVEST: {
+        // Harvest can both remove biomass and move biomass to the soil/litter
+        // pools
+        const HarvestParams *harvParams = event->eventParams;
+        const double fracRA = harvParams->fractionRemovedAbove;
+        const double fracRB = harvParams->fractionRemovedBelow;
+        const double fracTA = harvParams->fractionTransferredAbove;
+        const double fracTB = harvParams->fractionTransferredBelow;
+
+        // No need to allocate to biomass N pools, we don't track that N
+        // explicitly. We do need to handle soil and litter N, though.
+        // Note: ctx.nitrogenCycle implies ctx.litterPool
+        // Litter N increase
+        const double totalAbove = (envi.plantLeafC / params.leafCN) +
+                                  (envi.plantWoodC / params.woodCN);
+        const double totalBelow = (envi.fineRootC / params.fineRootCN) +
+                                  (envi.coarseRootC / params.woodCN);
+        double litterNAdd = fracTA * totalAbove;
+        double soilNAdd = fracTB * totalBelow;
+        fluxes.eventSoilOrgN += soilNAdd / climLen;
+        fluxes.eventLitterN += litterNAdd / climLen;
+
+        // MASS BALANCE: removed fractions are system outputs
+        // just plantWoodC here, not woodC
+        double outputN = (envi.plantWoodC / params.woodCN +
+                          envi.plantLeafC / params.leafCN) *
+                             fracRA +
+                         (envi.fineRootC / params.fineRootCN +
+                          envi.coarseRootC / params.woodCN) *
+                             fracRB;
+        fluxes.eventOutputN += outputN / climLen;
+
+        // clang-format off
+        appendLog(
+            event, 3,
+            "eventSoilOrgN", soilNAdd,
+            "eventLitterN", litterNAdd,
+            "eventOutputN", outputN);
+        // clang-format on
+      } break;
+      case TILLAGE: {
+        // Nothing to do for tillage events
+      } break;
+      case FERTILIZATION: {
+        const FertilizationParams *fertParams = event->eventParams;
+        double orgN = fertParams->orgN;
+        double minN = fertParams->minN;
+        // As the warning says in readEventData(), we ignore N when the
+        // nitrogen cycle model is off
+        // Implies ctx.litterPool
+        fluxes.eventLitterN += orgN / climLen;
+        fluxes.eventMinN += minN / climLen;
+
+        // MASS BALANCE: this is a system input
+        fluxes.eventInputN += (orgN + minN) / climLen;
+
+        // clang-format off
+        appendLog(event, 3,
+          "eventMinN", minN,
+          "eventLitterN", orgN,
+          "eventInputN", (orgN + minN));
+        // clang-format on
+      } break;
+      case LEAFON: {
         // Nitrogen is handled implicitly by relative CN ratios. Missing N
         // from low-N wood to higher-N leaves is accounted for in
         // calcNFixationAndUptakeFluxes() via calcPlantNDemandFlux()
 
         // 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
+        // eventInputN
       } break;
       case LEAFOFF: {
-        double leafOff = envi.plantLeafC * params.fracLeafFall;
-        fluxes.eventLeafOffLitter += leafOff / climLen;
-
-        double litterNAdd = 0.0;
-        double leafNResorption = 0.0;
-        if (ctx.nitrogenCycle) {
-          // Nitrogen - need to account for leaf N moving to litter, as with
-          // harvests
-          double leafN = leafOff / params.leafCN;
-          leafNResorption = leafN * params.leafNResorptionFrac;
-          litterNAdd = leafN - leafNResorption;
-          fluxes.eventLeafOffNResorption += leafNResorption / climLen;
-          fluxes.eventLitterN += litterNAdd / climLen;
+        // We need to use fluxes.eventLeafOffLitter here instead of
+        // recalculating it, as it may have been reduced. Note that this means
+        // we are processing ALL leaf-off events in one shot in the unlikely
+        // event that there are more than one in this time step.
+        if (!firstLeafOff) {
+          logInfo(
+              "Ignoring nitrogen effects of second (or more) leaf-off "
+              "event at year %d day %d; all nitrogen effects are recorded with "
+              "the first leaf-off event in this time step\n",
+              event->year, event->day);
+          break;
         }
+        firstLeafOff = 0;
+        // Nitrogen - need to account for leaf N moving to litter, as with
+        // harvests
+        double leafOff = fluxes.eventLeafOffLitterC * climLen;
+        double preResorp = fluxes.eventLeafOffNResorption;
+        double preLitter = fluxes.eventLeafOffLitterN;
+        calcLeafOffNEffects(leafOff, &fluxes.eventLeafOffNResorption,
+                            &fluxes.eventLeafOffLitterN);
+
+        double leafNResorptionFlux = fluxes.eventLeafOffNResorption - preResorp;
+        double litterNAddFlux = fluxes.eventLeafOffLitterN - preLitter;
 
         // clang-format off
-        writeEventOut(gEvent, 3,
-          "eventLeafOffLitter", leafOff,
-          "eventLeafOffNResorption", leafNResorption,
-          "eventLitterN", litterNAdd);
+        appendLog(event, 2,
+          "eventLeafOffNResorption", leafNResorptionFlux * climLen,
+          "eventLitterN", litterNAddFlux * climLen);
         // clang-format on
       } break;
       case PLANTDEATH:
-        // There should be no way to get here, but covering our bases...
-        logWarning("PLANTDEATH event found for year %d day %d, but not "
-                   "implemented as an input event; ignoring\n",
-                   gEvent->year, gEvent->day);
+        // Nothing to do here
         break;
       default:
-        logError("Unknown event type (%d) in processEvents()\n", gEvent->type);
+        logError("Unknown event type (%d) in processEvents()\n", event->type);
         exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM);
     }
 
-    gEvent = gEvent->nextEvent;
+    event = event->nextEvent;
   }
 }
 
@@ -745,6 +895,7 @@ void updatePoolsForEvents(void) {
   // CARBON
   // Harvest and planting events
   envi.plantWoodC += fluxes.eventWoodC * climate->length;
+  envi.plantCAccountingDelta += fluxes.eventAccountingC * climate->length;
   envi.plantLeafC += fluxes.eventLeafC * climate->length;
 
   // Harvest and fertilization events
@@ -759,12 +910,12 @@ void updatePoolsForEvents(void) {
   double eventLeafOnCreationFromRoot =
       fluxes.eventLeafOnCreation - fluxes.eventLeafOnCreationFromWood;
   envi.coarseRootC -= eventLeafOnCreationFromRoot * climate->length;
-  envi.plantLeafC += (fluxes.eventLeafOnCreation - fluxes.eventLeafOffLitter) *
+  envi.plantLeafC += (fluxes.eventLeafOnCreation - fluxes.eventLeafOffLitterC) *
                      climate->length;
   if (ctx.litterPool) {
-    envi.litterC += fluxes.eventLeafOffLitter * climate->length;
+    envi.litterC += fluxes.eventLeafOffLitterC * climate->length;
   } else {
-    envi.soilC += fluxes.eventLeafOffLitter * climate->length;
+    envi.soilC += fluxes.eventLeafOffLitterC * climate->length;
   }
 
   // Harvest and planting events
@@ -782,7 +933,8 @@ void updatePoolsForEvents(void) {
   if (ctx.nitrogenCycle) {
     envi.minN += fluxes.eventMinN * climate->length;
     envi.soilOrgN += fluxes.eventSoilOrgN * climate->length;
-    envi.litterN += fluxes.eventLitterN * climate->length;
+    envi.litterN +=
+        (fluxes.eventLitterN + fluxes.eventLeafOffLitterN) * climate->length;
     double leafOnNFlux = calcLeafOnNFromC(fluxes.eventLeafOnCreation);
     envi.plantStorageN +=
         (fluxes.eventLeafOffNResorption - leafOnNFlux) * climate->length;
@@ -799,6 +951,9 @@ void freeEventList(void) {
     if (prev->eventParams != NULL) {
       free(prev->eventParams);
     }
+    if (prev->logLine != NULL) {
+      dsFree(prev->logLine);
+    }
     free(prev);
   }
 }
diff --git a/src/sipnet/events.h b/src/sipnet/events.h
index 22a43fcb..ed3aa6fd 100644
--- a/src/sipnet/events.h
+++ b/src/sipnet/events.h
@@ -1,6 +1,8 @@
 #ifndef EVENTS_H
 #define EVENTS_H
 
+#include "common/util.h"
+
 typedef enum EventType {
   FERTILIZATION,
   HARVEST,
@@ -81,6 +83,8 @@ struct EventNode {
   event_type_t type;
   int year, day;
   void *eventParams;
+  int numLogParamPairs;
+  DynamicString *logLine;
   EventNode *nextEvent;
 };
 
@@ -117,29 +121,20 @@ EventNode *readEventData(const char *eventFile);
 void openEventOutFile(const char *eventOutFile, int printHeader);
 
 /*!
- * \brief Write a line to the event output file for a single oneEvent
- *
- * Writes a single oneEvent to the configured event output file. This is a
- * variadic function which expects to receive 2*numParams values in (char*,
- * double) pairs after the
- * numParams argument.
- *
- * Output format:
- *
- * year day event_type \<param_name>=\<delta>[,\<param_name>=\<delta>,...]
- *
- * \param oneEvent     Pointer to oneEvent node
- * \param numParams Number of param/value PAIRS to write
- * \param ...       Pairs of (char*, double) arguments to write, 2*numParams
- *                  values
+ * Append to an event's log line
  */
-void writeEventOut(EventNode *oneEvent, int numParams, ...);
+void appendLog(EventNode *event, int numParams, ...);
+
+/*!
+ * Write out all events for this time step
+ */
+void writeEventsOut(void);
 
 /*!
  * \brief Write a line to the event output file for a computed event
  *
- * Same as writeEventOut, but for events that are computed internally, such
- * as leaf on/leaf off events.
+ * Write an event that is computed internally, such as leaf on/leaf off or
+ * plant death events.
  *
  * Output format:
  *
@@ -187,17 +182,30 @@ void setupEvents(void);
 int isFirstEventBefore(int year, int day);
 
 /*!
- * \brief Process events for current location/year/day
+ * Return today's first event for SIPNET's multi-pass calls
+ */
+EventNode *getCurrentEvent(void);
+
+/*!
+ * \brief Process carbon effects from events for current day
  *
- * For a given year and day (as determined by the global `climate`
- * pointer), process all events listed in the global `events` pointer for the
- * referenced location.
+ * Process all events for the current day, calculating all carbon
+ * effects. For each event, modify flux variables according to the model for
+ * that event type.
+ */
+void processEventsForCarbon(EventNode *event);
+
+/*!
+ * \brief Process nitrogen effects from events for current day
  *
- * For each event, modify flux variables according to the model for that event,
- * and write a row to the configured event output file listing the modified
- * variables and the delta applied.
+ * Process all events for the current day, calculating all nitrogen
+ * effects. For each event, modify flux variables according to the model for
+ * that event type.
+ *
+ * Carbon and nitrogen effects are calculated separately to allow carbon
+ * limitation checks to be run before any nitrogen calculations are made.
  */
-void processEvents(void);
+void processEventsForNitrogen(EventNode *event);
 
 /*!
  * Update relevant environment pools after event fluxes have been calculated
@@ -209,6 +217,15 @@ void updatePoolsForEvents(void);
  */
 void freeEventList(void);
 
+typedef struct HarvestTrackersStruct {
+  double totalFracRemoved;
+  double totalFracTransferred;
+  double totalFracRemovedAbove;
+  double totalFracRemovedBelow;
+  double totalFracTransferredAbove;
+  double totalFracTransferredBelow;
+} HarvestTrackers;
+
 // Variables to track events with lingering effects
 typedef struct EventTrackerStruct {
   // Tillage effect on Rh; exponentially decays at each time step by a factor
@@ -216,8 +233,7 @@ typedef struct EventTrackerStruct {
   double d_till_mod;
 
   // Fraction removed and transferred from harvest event this time step
-  double harvestFracRemoved;
-  double harvestFracTransferred;
+  HarvestTrackers harvestTrackers;
 } EventTrackers;
 
 extern EventTrackers eventTrackers;
diff --git a/src/sipnet/limitations.c b/src/sipnet/limitations.c
index 726cedc2..6a5d3a9e 100644
--- a/src/sipnet/limitations.c
+++ b/src/sipnet/limitations.c
@@ -63,16 +63,87 @@ void checkLeafOnLimitation(double *leafOnFlux) {
   }
 }
 
+/**
+ * Check that negative growth is not driving a pool to end negative
+ *
+ * Adjust if necessary
+ */
+static void checkNegativeCreation(void) {
+  // In the case of negative growth (mean npp < 0), we might be allocating that
+  // negative growth to a pool that can't handle it (e.g., leaf creation is
+  // negative, but leaf pool is already at 0). In those cases, adjust
+  // appropriately.
+
+  double len = climate->length;
+
+  // Above ground
+  // If leafCreation is too negative, we need to deduct from wood instead
+  // Need to make sure we don't go too far the other way. Leaf off litter
+  // (either fluxes.leafLitter or fluxes.eventLeafOffLitter) might also
+  // incorrectly drive the pool negative, sp adjust for that too.
+
+  // First we handle leaf litter. Assuming params.leafTurnoverRate is valid
+  // (ie, <=1), availableLeafRate should be non-negative.
+  double availableLeafRate = envi.plantLeafC / len - fluxes.leafLitter;
+
+  // Reminder: litter and creation fluxes "point" in the opposite direction, so
+  // their signs are opposite here
+  double deficit = fluxes.leafOffLitter + fluxes.eventLeafOffLitterC -
+                   fluxes.leafCreation - availableLeafRate;
+  if (deficit > 0) {
+    // If negative growth + leaf off is too much, let's first reduce leaf off
+    if (fluxes.eventLeafOffLitterC > 0) {
+      double adjust = fmin(fluxes.eventLeafOffLitterC, deficit);
+      fluxes.eventLeafOffLitterC -= adjust;
+      deficit -= adjust;
+    }
+    if (deficit > 0) {
+      if (fluxes.leafOffLitter > 0) {
+        double adjust = fmin(fluxes.leafOffLitter, deficit);
+        fluxes.leafOffLitter -= adjust;
+        deficit -= adjust;
+      }
+    }
+    // Next, move some negative growth to the wood pool
+    if (deficit > 0) {
+      // If this is too much for the wood pool to handle, we have an error that
+      // will be caught in ensureNonNegative
+      fluxes.woodCreation -= deficit;
+      fluxes.leafCreation += deficit;
+    }
+  }
+
+  // Below ground
+  double fineRootDeficit =
+      envi.fineRootC / len + fluxes.fineRootCreation - fluxes.fineRootLoss;
+  double coarseRootDeficit = envi.coarseRootC / len +
+                             fluxes.coarseRootCreation - fluxes.coarseRootLoss;
+  if ((fineRootDeficit < 0.0) != (coarseRootDeficit < 0.0)) {
+    // If neither are negative, nothing to do
+    // If both are negative, the plant will die in checkForMortality()
+    if (fineRootDeficit < 0.0) {
+      fluxes.coarseRootCreation += fineRootDeficit;
+      fluxes.fineRootCreation -= fineRootDeficit;
+    }
+    if (coarseRootDeficit < 0.0) {
+      fluxes.fineRootCreation += coarseRootDeficit;
+      fluxes.coarseRootCreation -= coarseRootDeficit;
+    }
+  }
+}
+
 /**
  * Check for nitrogen limitation, and reduce growth if needed
  */
 static void checkNitrogenLimitation(void) {
+  double len = climate->length;
+
   // First, determine if we are in a nitrogen-limited situation. The uptake
   // flux has already taken the storage pool into account, so we just need to
   // see if that uptake is too much, taking into account other fluxes to the
   // minN pool.
   // Calc total delta to minN pool
-  double len = climate->length;
+  // double len = climate->length;
   double uptakeDemand = fluxes.nUptake * len;
   double nonUptakeDelta = calcMinNNonUptakeFluxes() * len;
   double availableMinN = envi.minN + nonUptakeDelta;
@@ -108,8 +179,22 @@ static void checkNitrogenLimitation(void) {
     fluxes.fineRootCreation *= reduction;
     fluxes.coarseRootCreation *= reduction;
 
-    // Reset fixation and uptake
-    calcNFixationAndUptakeFluxes();
+    // If there was a leaf-off event in this step, it is possible that we have
+    // reduced leaf growth down to the point where the pool will go negative
+    // Re-call that check
+    checkNegativeCreation();
+
+    // Reset and recalc event leaf-off N
+    // Nitrogen for "calculated" leaf-off events will be redone in
+    // calcNitrogenFluxes below
+    fluxes.eventLeafOffNResorption = 0.0;
+    fluxes.eventLeafOffLitterN = 0.0;
+    calcLeafOffNEffects(fluxes.eventLeafOffLitterC * len,
+                        &fluxes.eventLeafOffNResorption,
+                        &fluxes.eventLeafOffLitterN);
+
+    // Reset N calculations
+    calcNitrogenFluxes();
   }
 }
 
@@ -138,48 +223,5 @@ void checkLimitations(void) {
   }
 }
 
-/**
- * Check that negative growth is not driving a pool to end negative
- *
- * Adjust if necessary
- */
-static void checkNegativeCreation(void) {
-  // In the case of negative growth (mean npp < 0), we might be allocating that
-  // negative growth to a pool that can't handle it (e.g., leaf creation is
-  // negative, but leaf pool is already at 0). In those cases, adjust
-  // appropriately.
-
-  double len = climate->length;
-  // Above ground
-  // If leafCreation is too negative, we need to deduct from wood instead
-  // Use only the continuous turnover term to match previous logic - but see
-  // SIPNET issue #372.
-  double leafLitterTurnover = envi.plantLeafC * params.leafTurnoverRate;
-  double leafDeficit =
-      envi.plantLeafC / len + fluxes.leafCreation - leafLitterTurnover;
-  if (leafDeficit < 0) {
-    fluxes.woodCreation += leafDeficit;
-    fluxes.leafCreation -= leafDeficit;
-  }
-
-  // Below ground
-  double fineRootDeficit =
-      envi.fineRootC / len + fluxes.fineRootCreation - fluxes.fineRootLoss;
-  double coarseRootDeficit = envi.coarseRootC / len +
-                             fluxes.coarseRootCreation - fluxes.coarseRootLoss;
-  if ((fineRootDeficit < 0.0) != (coarseRootDeficit < 0.0)) {
-    // If neither are negative, nothing to do
-    // If both are negative, the plant will die in checkForMortality()
-    if (fineRootDeficit < 0.0) {
-      fluxes.coarseRootCreation += fineRootDeficit;
-      fluxes.fineRootCreation -= fineRootDeficit;
-    }
-    if (coarseRootDeficit < 0.0) {
-      fluxes.fineRootCreation += coarseRootDeficit;
-      fluxes.coarseRootCreation -= coarseRootDeficit;
-    }
-  }
-}
-
 // See limitations.h
 void checkCarbonLimitations(void) { checkNegativeCreation(); }
diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c
index 739eb290..70ca096e 100644
--- a/src/sipnet/nitrogen.c
+++ b/src/sipnet/nitrogen.c
@@ -68,7 +68,7 @@ static void calcNPoolFluxes(void) {
   // litter (modified by leaf N resorption), and N loss due to mineralization.
   // N added via fertilization is handled elsewhere.
   fluxes.nOrgLitter =
-      fluxes.leafLitter / params.leafCN - fluxes.leafOffNResorption +
+      getLeafLitterFlux() / params.leafCN - fluxes.leafOffNResorption +
       fluxes.woodLitter / params.woodCN - litterMin -
       fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction);
 
@@ -186,12 +186,11 @@ void calcNResorptionFluxes(void) {
          fluxes.fineRootCreation / params.fineRootCN);
   }
 
-  // Leaf litter resorption; at this point, fluxes.leafLitter counts both normal
-  // turnover and leaf-off calcs. Note that event leaf off is handled in
-  // events.c
+  // Leaf litter resorption; at this point. Note that event leaf off is handled
+  // in events.c
   double nResorp =
-      params.leafNResorptionFrac * fluxes.leafLitter / params.leafCN;
-  fluxes.leafOffNResorption += nResorp;
+      params.leafNResorptionFrac * getLeafLitterFlux() / params.leafCN;
+  fluxes.leafOffNResorption = nResorp;
 
   // TODO: Should we resorb N from wood litter?
 }
@@ -238,3 +237,13 @@ void updateNitrogenPools(void) {
   // Litter organic N
   envi.litterN += fluxes.nOrgLitter * climate->length;
 }
+
+void calcLeafOffNEffects(double leafOffC, double *resorptionFlux,
+                         double *litterFlux) {
+  double climLen = climate->length;
+  double leafN = leafOffC / params.leafCN;
+  double leafNResorption = leafN * params.leafNResorptionFrac;
+  double litterNAdd = leafN - leafNResorption;
+  *resorptionFlux += leafNResorption / climLen;
+  *litterFlux += litterNAdd / climLen;
+}
diff --git a/src/sipnet/nitrogen.h b/src/sipnet/nitrogen.h
index db7a475b..8a30830a 100644
--- a/src/sipnet/nitrogen.h
+++ b/src/sipnet/nitrogen.h
@@ -78,4 +78,13 @@ void updateNitrogenPools(void);
  */
 void calcNResorptionFluxes(void);
 
+/**
+ * Helper function for handling leaf-off events
+ *
+ * This is used by both processEventsForNitrogen and the nitrogen limitation
+ * check in limitations.
+ */
+void calcLeafOffNEffects(double leafOffC, double *resorptionFlux,
+                         double *litterFlux);
+
 #endif  // NITROGEN_H
diff --git a/src/sipnet/restart.c b/src/sipnet/restart.c
index e8c2a4b6..78689e22 100644
--- a/src/sipnet/restart.c
+++ b/src/sipnet/restart.c
@@ -37,7 +37,7 @@
 #define NUM_TRACKER_FIELDS 33
 #define NUM_PHENOLOGY_TRACKERS_FIELDS 3
 #define NUM_SURVIVAL_TRACKERS_FIELDS 1
-#define NUM_EVENT_TRACKERS_FIELDS 3
+#define NUM_EVENT_TRACKERS_FIELDS 7
 
 // This one shouldn't change
 #define NUM_END_FIELDS 1
@@ -291,9 +291,13 @@ void initResetState(RestartState *state, MeanTracker *npp) {
   }
 
   ind = 0;
-  state->eventPF[ind++] = (StateField){"event_trackers.d_till_mod",             FT_DOUBLE,  &eventTrackers.d_till_mod,             0};
-  state->eventPF[ind++] = (StateField){"event_trackers.harvestFracRemoved",     FT_DOUBLE,  &eventTrackers.harvestFracRemoved,     0};
-  state->eventPF[ind++] = (StateField){"event_trackers.harvestFracTransferred", FT_DOUBLE,  &eventTrackers.harvestFracTransferred, 0};
+  state->eventPF[ind++] = (StateField){"event_trackers.d_till_mod",                                FT_DOUBLE,  &eventTrackers.d_till_mod,                                0};
+  state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracRemoved",          FT_DOUBLE,  &eventTrackers.harvestTrackers.totalFracRemoved,          0};
+  state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracTransferred",      FT_DOUBLE,  &eventTrackers.harvestTrackers.totalFracTransferred,      0};
+  state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracRemovedAbove",     FT_DOUBLE,  &eventTrackers.harvestTrackers.totalFracRemovedAbove,     0};
+  state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracRemovedBelow",     FT_DOUBLE,  &eventTrackers.harvestTrackers.totalFracRemovedBelow,     0};
+  state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracTransferredAbove", FT_DOUBLE,  &eventTrackers.harvestTrackers.totalFracTransferredAbove, 0};
+  state->eventPF[ind++] = (StateField){"event_trackers.harvestTrackers.totalFracTransferredBelow", FT_DOUBLE,  &eventTrackers.harvestTrackers.totalFracTransferredBelow, 0};
   state->eventPF[ind++] = (StateField){"event_trackers.invalid",                FT_INVALID, NULL, FIELD_INVALID};
   if (ind != NUM_EVENT_TRACKERS_FIELDS + 1) {
     logInternalError("Restart array size mismatch: eventPF\n");
diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c
index 5b3a86b7..53de57de 100644
--- a/src/sipnet/sipnet.c
+++ b/src/sipnet/sipnet.c
@@ -791,19 +791,11 @@ void calcWoodAndLeafFluxes(void) {
  * the start or end of the growing season. These transition fluxes are tracked
  * separately from the continuous leaf creation and litter fluxes calculated in
  * calcLeafFluxes().
- *
- * @param[out] leafOnCreation Additional leaf creation flux at leaf-on
- *   (g C/m^2 ground/day)
- * @param[out] leafOnFromWood Carbon transferred from wood to leaves at
- *   leaf-on (g C/m^2 ground/day)
- * @param[out] leafLitter Additional leaf litter flux at leaf-off
- *   (g C/m^2 ground/day)
- * @param[in] plantLeafC Leaf carbon pool size (g C/m^2 ground area)
  */
-void calcLeafOnOffFluxes(double *leafOnCreation, double *leafOnFromWood,
-                         double *leafLitter, double plantLeafC) {
+void calcLeafOnOffFluxes(void) {
   // Calc additional fluxes at start/end of growing season
   // Note that these are basically events, and we will track them as such
+  double len = climate->length;
 
   // first check for new year; if new year, reset trackers (since we haven't
   // done leaf growth or fall yet in this new year):
@@ -816,25 +808,30 @@ void calcLeafOnOffFluxes(double *leafOnCreation, double *leafOnFromWood,
   // check for start of growing season:
   if (!phenologyTrackers.didLeafGrowth && pastLeafGrowth()) {
     // we just reached the start of the growing season
-    double leafOn = params.leafGrowth / climate->length;
+    double leafOn = params.leafGrowth / len;
+    double leafOnFromWood = 0;
     checkLeafOnLimitation(&leafOn);
-    *leafOnCreation += leafOn;
+    fluxes.leafOnCreation += leafOn;
     double totalSourceC = envi.plantWoodC + envi.coarseRootC;
     if (totalSourceC > TINY) {
-      *leafOnFromWood += leafOn * envi.plantWoodC / totalSourceC;
+      leafOnFromWood = leafOn * envi.plantWoodC / totalSourceC;
+      fluxes.leafOnCreationFromWood += leafOnFromWood;
     }
     phenologyTrackers.didLeafGrowth = 1;
-    // This is a computed event - however, the value may get reduced by
-    // nitrogen limitation. The writeEvent call is in writeLeafOnEventIfNeeded,
-    // called after N limiting is checked.
+
+    if (ctx.events && leafOn > TINY) {
+      writeComputedEventOut(climate->year, climate->day,
+                            eventTypeToString(LEAFON), 2, "leafOnCreation",
+                            leafOn * len, "leafOnCreationFromWood",
+                            leafOnFromWood * len);
+    }
   }
 
   // check for end of growing season:
   if (!phenologyTrackers.didLeafFall && pastLeafFall()) {
     // we just reached the end of the growing season
-    double len = climate->length;
-    double leafOff = (plantLeafC * params.fracLeafFall) / len;
-    *leafLitter += leafOff;
+    double leafOff = (envi.plantLeafC * params.fracLeafFall) / len;
+    fluxes.leafOffLitter += leafOff;
     phenologyTrackers.didLeafFall = 1;
     if (leafOff > TINY && ctx.events) {
       writeComputedEventOut(climate->year, climate->day,
@@ -1224,31 +1221,6 @@ void calcMethaneFlux(void) {
  */
 void resetFluxes(void) { fluxes = (struct FluxVars){0}; }
 
-/**
- * Write out a leaf-on event if one happened
- *
- * Delayed event writing for leaf-on, if appropriate, since the value may
- * have changed due to N limitation
- */
-void writeLeafOnEventIfNeeded(void) {
-  const char *type = eventTypeToString(LEAFON);
-  const double len = climate->length;
-  if (fluxes.leafOnCreation > TINY && ctx.events) {
-    writeComputedEventOut(climate->year, climate->day, type, 2,
-                          "leafOnCreation", fluxes.leafOnCreation * len,
-                          "leafOnCreationFromWood",
-                          fluxes.leafOnCreationFromWood * len);
-  }
-  if (fluxes.eventLeafOnCreation > TINY && ctx.events) {
-    // Not really a computed event, but we don't have the event object here, so
-    // we use this mechanism
-    writeComputedEventOut(
-        climate->year, climate->day, type, 2, "eventLeafOnCreation",
-        fluxes.eventLeafOnCreation * len, "eventLeafOnCreationFromWood",
-        fluxes.eventLeafOnCreationFromWood * len);
-  }
-}
-
 /*!
  * Calculate flux terms for sipnet as part of main model flow
  *
@@ -1276,6 +1248,9 @@ void calculateFluxes(void) {
   // Psn, moisture and water fluxes
   lai = envi.plantLeafC / params.leafCSpWt;  // current lai
 
+  // Today's first event
+  EventNode *event = getCurrentEvent();
+
   potPsn(&potGrossPsn, &baseFolResp, lai, climate->tair, climate->vpd,
          climate->par);
   moisture(&(fluxes.transpiration), &dWater, potGrossPsn, climate->vpd,
@@ -1288,6 +1263,9 @@ void calculateFluxes(void) {
                       fluxes.snowMelt, fluxes.transpiration);
   getGpp(&(fluxes.photosynthesis), potGrossPsn, dWater);
 
+  // First pass for events: carbon effects
+  processEventsForCarbon(event);
+
   // Vegetation respiration
   if (ctx.growthResp) {
     vegResp2(&folResp, &woodResp, &growthResp, baseFolResp);
@@ -1301,8 +1279,7 @@ void calculateFluxes(void) {
   calcWoodAndLeafFluxes();
 
   // Leaf on/off
-  calcLeafOnOffFluxes(&fluxes.leafOnCreation, &fluxes.leafOnCreationFromWood,
-                      &fluxes.leafLitter, envi.plantLeafC);
+  calcLeafOnOffFluxes();
 
   // Litter pool, if LITTER is on
   calcLitterFluxes();
@@ -1327,15 +1304,17 @@ void calculateFluxes(void) {
   // so make sure this stays at the bottom of this function (or after the
   // carbon and water calcs, at least).
   if (ctx.nitrogenCycle) {
+    // Second pass for events: nitrogen effects
+    processEventsForNitrogen(event);
+    // General nitrogen fluxes
     calcNitrogenFluxes();
   }
 
   // Check limitations; must be after all flux calcs
   checkLimitations();
 
-  // Delayed write needed due to N limitation check - leafOn value may have
-  // changed
-  writeLeafOnEventIfNeeded();
+  // Write to events file for all of today's events
+  writeEventsOut();
 }
 
 // /////////////////////// //
@@ -1478,7 +1457,7 @@ void updateTrackers(double oldSoilWater) {
 
   // If we get another event flux in this function, we should create an
   // updateTrackersForEvents() function in events.c|h
-  trackers.yearlyLitter += fluxes.leafLitter + fluxes.eventLeafOffLitter;
+  trackers.yearlyLitter += getLeafLitterFlux() + fluxes.eventLeafOffLitterC;
 
   if (ctx.gdd) {
     trackers.gdd += climate->gdd;
@@ -1529,8 +1508,16 @@ void initPhenologyTrackers(void) {
                                                // this year
 }
 
-// Check that woodC and total root C are both positive
+// Check that woodC and total root C are both positive and that there was no
+// terminating harvest
 int hasSufficientBiomass(void) {
+  // If there was a harvest termination event, the answer is no
+  if (eventTrackers.harvestTrackers.totalFracRemoved +
+          eventTrackers.harvestTrackers.totalFracTransferred >
+      1.0 - TINY) {
+    return 0;
+  }
+
   double totalWoodC = getTotalWoodC();
   double totalRootC = envi.fineRootC + envi.coarseRootC;
   // We want to check that both plantWoodC AND totalWoodC are positive, as well
@@ -1608,7 +1595,7 @@ void updateMainPools(void) {
   //     L_L = fluxes.leafLitter
   // Note: we have split leafCreation into two parts
   envi.plantLeafC +=
-      (fluxes.leafCreation + fluxes.leafOnCreation - fluxes.leafLitter) *
+      (fluxes.leafCreation + fluxes.leafOnCreation - getLeafLitterFlux()) *
       climate->length;
 
   // :: from [1], eq (A4), where:
@@ -1649,7 +1636,7 @@ void updatePoolsForSoil(void) {
         ctx.carbonSaturation ? unitClip(envi.soilC / params.soilCSaturation)
                              : 0.0;
     // :: from [2], litter model description
-    envi.litterC += (fluxes.woodLitter + fluxes.leafLitter +
+    envi.litterC += (fluxes.woodLitter + getLeafLitterFlux() +
                      (soilInputs * saturationFraction) - fluxes.litterToSoil -
                      fluxes.rLitter - fluxes.litterMethane) *
                     climate->length;
@@ -1661,12 +1648,12 @@ void updatePoolsForSoil(void) {
     // Normal pool (single pool, no microbes)
     // :: from [1] (and others, TBD), eq (A3), where:
     //     L_w = fluxes.woodLitter
-    //     L_l = fluxes.leafLitter
+    //     L_l = fluxes.leafLitter + fluxes.leafOffLitter
     //     R_h = fluxes.rSoil
     // :: from [3], root terms
     envi.soilC +=
         (fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.woodLitter +
-         fluxes.leafLitter - fluxes.rSoil - fluxes.soilMethane) *
+         getLeafLitterFlux() - fluxes.rSoil - fluxes.soilMethane) *
         climate->length;
   }
 
@@ -1706,18 +1693,15 @@ void checkForMortality(void) {
     plantSurvivalTracker.isAlive = 0;
     double totalWoodC = getTotalWoodC();
     double totalRootC = envi.fineRootC + envi.coarseRootC;
-
-    if (eventTrackers.harvestFracRemoved +
-            eventTrackers.harvestFracTransferred >=
-        TINY) {
+    HarvestTrackers *ht = &eventTrackers.harvestTrackers;
+    int harvestOccurred =
+        ht->totalFracRemoved + ht->totalFracTransferred >= TINY;
+    if (harvestOccurred) {
       logInfo("Plant mortality detected after harvest event: total fraction "
-              "removed %.3f total fraction transferred %.3f; woodC %f "
-              "totalWoodC %f coarseRootC %f fineRootC %f year %d day %d "
+              "removed %.3f total fraction transferred %.3f on year %d day %d "
               "time %6.3f; zeroing out biomass pools\n",
-              eventTrackers.harvestFracRemoved,
-              eventTrackers.harvestFracTransferred, envi.plantWoodC, totalWoodC,
-              envi.coarseRootC, envi.fineRootC, climate->year, climate->day,
-              climate->time);
+              ht->totalFracRemoved, ht->totalFracTransferred, climate->year,
+              climate->day, climate->time);
     } else {
       logWarning(
           "Plant mortality detected as wood or total root carbon is zero "
@@ -1732,19 +1716,34 @@ void checkForMortality(void) {
     // multiple processes being modeled, it's believable that there may be a bit
     // of overshoot when a plant dies - for example, a 100% harvest event with
     // any overall loss (respiration, turnover).
-    envi.soilC += totalRootC;
+
+    // Also, if there was a harvest, reduce by the appropriate removal fraction
+    // (For no harvest, these will be 1)
+    double aboveRemovalReduction = (1 - ht->totalFracRemovedAbove);
+    double belowRemovalReduction = (1 - ht->totalFracRemovedBelow);
+
+    envi.soilC += totalRootC * belowRemovalReduction;
+    double aboveC =
+        envi.plantWoodC + envi.plantLeafC + envi.plantCAccountingDelta;
     if (ctx.litterPool) {
-      envi.litterC +=
-          envi.plantWoodC + envi.plantLeafC + envi.plantCAccountingDelta;
+      envi.litterC += aboveC * aboveRemovalReduction;
     } else {
-      envi.soilC +=
-          envi.plantWoodC + envi.plantLeafC + envi.plantCAccountingDelta;
+      envi.soilC += aboveC * aboveRemovalReduction;
     }
+    fluxes.eventOutputC += (aboveC * ht->totalFracRemovedAbove +
+                            totalRootC * ht->totalFracRemovedBelow) /
+                           climate->length;
+
     if (ctx.nitrogenCycle) {  // litter pool implied
-      envi.soilOrgN +=
+      double aboveN =
+          envi.plantWoodC / params.woodCN + envi.plantLeafC / params.leafCN;
+      double belowN =
           envi.fineRootC / params.fineRootCN + envi.coarseRootC / params.woodCN;
-      envi.litterN += envi.plantWoodC / params.woodCN +
-                      envi.plantLeafC / params.leafCN + envi.plantStorageN;
+      envi.soilOrgN += belowN * belowRemovalReduction;
+      envi.litterN += aboveN * aboveRemovalReduction + envi.plantStorageN;
+      fluxes.eventOutputN += (aboveN * ht->totalFracRemovedAbove +
+                              belowN * ht->totalFracRemovedBelow) /
+                             climate->length;
     }
 
     // Force pools to zero
@@ -1760,11 +1759,11 @@ void checkForMortality(void) {
     resetMeanTracker(meanNPP, 0.0);
 
     if (ctx.events) {
-      writeComputedEventOut(
-          climate->year, climate->day, eventTypeToString(PLANTDEATH), 4,
-          "harvestFracRemoved", eventTrackers.harvestFracRemoved,
-          "harvestFracTransferred", eventTrackers.harvestFracTransferred,
-          "totalWoodC", totalWoodC, "totalRootC", totalRootC);
+      writeComputedEventOut(climate->year, climate->day,
+                            eventTypeToString(PLANTDEATH), 4,
+                            "harvestFracRemoved", ht->totalFracRemoved,
+                            "harvestFracTransferred", ht->totalFracTransferred,
+                            "totalWoodC", totalWoodC, "totalRootC", totalRootC);
     }
   }
 }
@@ -1791,12 +1790,13 @@ void updatePoolsAndBalance() {
     updateNitrogenPools();
   }
 
+  // Check for obvious plant death (wood and/or roots at zero); also
+  // correct for harvest termination
+  checkForMortality();
+
   // Calc total C and N after pool updates
   updateBalanceTrackerPostUpdate();
 
-  // Check for obvious plant death (wood and/or roots at zero)
-  checkForMortality();
-
   // Verify none of our stocks have gone negative (set any that are to zero).
   ensureNonNegativeStocks();
 
@@ -1833,12 +1833,7 @@ void updateState(void) {
   ///////////////////////
   // 1. Calculate Fluxes
 
-  // All event handling, which is modeled as fluxes. Note that we have this
-  // before the other fluxes so that everything is in place when we consider
-  // N limitation at the end of calculateFluxes().
-  processEvents();
-
-  // All non-event fluxes
+  // All modeled fluxes, including events
   calculateFluxes();
 
   ///////////////////////
diff --git a/src/sipnet/state.c b/src/sipnet/state.c
index fdaa4391..b082cbcd 100644
--- a/src/sipnet/state.c
+++ b/src/sipnet/state.c
@@ -17,3 +17,7 @@ PlantSurvivalTracker plantSurvivalTracker;
 double getTotalWoodC(void) {
   return envi.plantWoodC + envi.plantCAccountingDelta;
 }
+
+double getLeafLitterFlux(void) {
+  return fluxes.leafLitter + fluxes.leafOffLitter;
+}
diff --git a/src/sipnet/state.h b/src/sipnet/state.h
index 658e106f..73ef466f 100644
--- a/src/sipnet/state.h
+++ b/src/sipnet/state.h
@@ -473,7 +473,7 @@ typedef struct FluxVars {
 
   // GROSS photosynthesis (g C * m^-2 ground area * day^-1)
   double photosynthesis;
-  // Leaf fall (g C * m^-2 ground area * day^-1)
+  // Leaf fall (g C * m^-2 ground area * day^-1); excludes leaf-off events
   double leafLitter;
   // Wood flux to litter (g C * m^-2 ground area * day^-1)
   double woodLitter;
@@ -550,6 +550,9 @@ typedef struct FluxVars {
   // Portion of leaf-on creation C that comes from wood C (the rest comes from
   // coarse root C)
   double leafOnCreationFromWood;
+  // leaf fall from calculated leaf off 'events' (gdd, soil temp, day of year)
+  // (g C * m^-2 grount area * day^-1)
+  double leafOffLitter;
 
   // ****************************************
   // Fluxes for nitrogen cycle
@@ -595,6 +598,8 @@ typedef struct FluxVars {
   double eventLeafC;
   // plantWoodC addition
   double eventWoodC;
+  // N-free accounting carbon removed by harvest
+  double eventAccountingC;
   // plantFineRootC addition
   double eventFineRootC;
   // plantCoarseRootC addition
@@ -631,7 +636,9 @@ typedef struct FluxVars {
   // coarse root C)
   double eventLeafOnCreationFromWood;
   // Transfer from leafC to soil/litter from a leaf-off event
-  double eventLeafOffLitter;
+  double eventLeafOffLitterC;
+  // Transfer from leafN to soil/litter from a leaf-off event
+  double eventLeafOffLitterN;
   // Resorption of leaf N from a leaf-off event
   double eventLeafOffNResorption;
 
@@ -755,4 +762,6 @@ extern PlantSurvivalTracker plantSurvivalTracker;
 
 double getTotalWoodC(void);
 
+double getLeafLitterFlux(void);
+
 #endif  // SIPNET_STATE_H
diff --git a/tests/sipnet/test_events_infrastructure/events_output_header.out b/tests/sipnet/test_events_infrastructure/events_output_header.out
index ec34fc93..0bda11a3 100644
--- a/tests/sipnet/test_events_infrastructure/events_output_header.out
+++ b/tests/sipnet/test_events_infrastructure/events_output_header.out
@@ -1,7 +1,7 @@
 year  day  type     param_name=delta[,param_name=delta,...]
-2023   65  plant    eventLeafC=3.00,eventWoodC=4.00,eventFineRootC=5.00,eventCoarseRootC=6.00,eventInputC=18.00,eventInputN=0.00
+2023   65  plant    eventLeafC=3.00,eventWoodC=4.00,eventFineRootC=5.00,eventCoarseRootC=6.00,eventInputC=18.00
 2023   70  irrig    eventSoilWater=5.00,eventEvap=0.00
-2023  200  harv     eventSoilC=1.90,eventLitterC=3.56,eventLeafC=-5.93,eventWoodC=-4.75,eventFineRootC=-3.73,eventCoarseRootC=-3.89,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=12.83,eventOutputN=0.00
-2024   65  plant    eventLeafC=3.00,eventWoodC=5.00,eventFineRootC=7.00,eventCoarseRootC=9.00,eventInputC=24.00,eventInputN=0.00
+2023  200  harv     eventSoilC=1.90,eventLitterC=3.56,eventLeafC=-5.93,eventWoodC=-4.75,eventAccountingC=-0.00,eventFineRootC=-3.73,eventCoarseRootC=-3.89,eventOutputC=12.83
+2024   65  plant    eventLeafC=3.00,eventWoodC=5.00,eventFineRootC=7.00,eventCoarseRootC=9.00,eventInputC=24.00
 2024   70  irrig    eventSoilWater=2.50,eventEvap=2.50
-2024  200  harv     eventSoilC=2.74,eventLitterC=1.51,eventLeafC=-1.39,eventWoodC=-1.63,eventFineRootC=-2.52,eventCoarseRootC=-2.97,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=4.25,eventOutputN=0.00
+2024  200  harv     eventSoilC=2.74,eventLitterC=1.51,eventLeafC=-1.39,eventWoodC=-1.63,eventAccountingC=-0.00,eventFineRootC=-2.52,eventCoarseRootC=-2.97,eventOutputC=4.25
diff --git a/tests/sipnet/test_events_infrastructure/events_output_no_header.out b/tests/sipnet/test_events_infrastructure/events_output_no_header.out
index e8465063..553d7073 100644
--- a/tests/sipnet/test_events_infrastructure/events_output_no_header.out
+++ b/tests/sipnet/test_events_infrastructure/events_output_no_header.out
@@ -1,6 +1,6 @@
-2023   65  plant    eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00,eventInputN=0.00
+2023   65  plant    eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00
 2023   70  irrig    eventSoilWater=5.00,eventEvap=0.00
-2023  200  harv     eventSoilC=12.40,eventLitterC=0.00,eventLeafC=-4.80,eventWoodC=-3.20,eventFineRootC=-4.80,eventCoarseRootC=-4.80,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=5.20,eventOutputN=0.00
-2024   65  plant    eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00,eventInputN=0.00
+2023  200  harv     eventSoilC=12.40,eventLitterC=0.00,eventLeafC=-4.80,eventWoodC=-3.20,eventAccountingC=-0.00,eventFineRootC=-4.80,eventCoarseRootC=-4.80,eventOutputC=5.20
+2024   65  plant    eventLeafC=10.00,eventWoodC=5.00,eventFineRootC=4.00,eventCoarseRootC=3.00,eventInputC=22.00
 2024   70  irrig    eventSoilWater=2.50,eventEvap=2.50
-2024  200  harv     eventSoilC=12.14,eventLitterC=0.00,eventLeafC=-10.32,eventWoodC=-5.88,eventFineRootC=-2.88,eventCoarseRootC=-2.48,eventSoilOrgN=0.00,eventLitterN=0.00,eventOutputC=9.42,eventOutputN=0.00
+2024  200  harv     eventSoilC=12.14,eventLitterC=0.00,eventLeafC=-10.32,eventWoodC=-5.88,eventAccountingC=-0.00,eventFineRootC=-2.88,eventCoarseRootC=-2.48,eventOutputC=9.42
diff --git a/tests/sipnet/test_events_types/testEventHarvest.c b/tests/sipnet/test_events_types/testEventHarvest.c
index 9dbc98f4..1cb089fd 100644
--- a/tests/sipnet/test_events_types/testEventHarvest.c
+++ b/tests/sipnet/test_events_types/testEventHarvest.c
@@ -140,9 +140,39 @@ int run(void) {
   return status;
 }
 
+int checkAccountingHarvest(void) {
+  int status = 0;
+  for (int i = -1; i <= 1; i++) {
+    prepTypesTest();
+    updateIntContext("litterPool", 1, CTX_TEST);
+    updateIntContext("nitrogenCycle", 1, CTX_TEST);
+    initEnv();
+    envi.plantCAccountingDelta = i;
+    updateBalanceTrackerPreUpdate();
+    initEvents("events_one_harvest.in", "events.out", 0);
+    setupEvents();
+    procEvents();
+    closeEventOutFile();
+    updateBalanceTrackerPostUpdate();
+    const double dc = balanceTracker.postTotalC - balanceTracker.preTotalC +
+                      fluxes.eventOutputC * climate->length;
+    const double dn = balanceTracker.postTotalN - balanceTracker.preTotalN +
+                      fluxes.eventOutputN * climate->length;
+    if (fabs(dc) > 1e-12 || fabs(dn) > 1e-12 ||
+        fabs(envi.plantWoodC - 1.8) > 1e-12 ||
+        fabs(envi.plantCAccountingDelta - 0.6 * i) > 1e-12) {
+      logTest("Accounting harvest failed: A=%d, C residual=%.15g, N "
+              "residual=%.15g\n",
+              i, dc, dn);
+      status = 1;
+    }
+  }
+  return status;
+}
+
 int main(void) {
   logTest("Starting run()\n");
-  int status = run();
+  int status = run() | checkAccountingHarvest();
   if (status) {
     logTest("FAILED testEventHarvest with status %d\n", status);
     exit(status);
diff --git a/tests/sipnet/test_modeling/Makefile b/tests/sipnet/test_modeling/Makefile
index 5cb4a25b..5ab3bfb5 100644
--- a/tests/sipnet/test_modeling/Makefile
+++ b/tests/sipnet/test_modeling/Makefile
@@ -8,7 +8,16 @@ LDFLAGS=-L$(ROOT_DIR)/libs
 LDLIBS=-lsipnet -lsipnet_common -lm
 
 # List test files in this directory here
-TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c testCarbonSaturation.c testPlantMortality.c testFluxCalculations.c
+TEST_CFILES= \
+  testBalance.c \
+  testCarbonSaturation.c \
+  testCompleteHarvest.c \
+  testDependencyFunctions.c \
+  testFluxCalculations.c \
+  testMethane.c \
+  testNitrogenCycle.c \
+  testPlantMortality.c \
+  testSoilMoisture.c \
 
 # The rest is boilerplate, likely copyable as is to a new test directory
 TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o)
diff --git a/tests/sipnet/test_modeling/termination.clim b/tests/sipnet/test_modeling/termination.clim
new file mode 100644
index 00000000..79baf35a
--- /dev/null
+++ b/tests/sipnet/test_modeling/termination.clim
@@ -0,0 +1,3 @@
+2017	19	21	0.125	10.574335	8.893	6.502	3.059	256.369	119.4744	1025.2	3.23371
+2017	20	0	0.125	10.85354	8.879	5.898	2.163	281.998	119.9049	1023.9	2.411573
+2017	20	3	0.125	8.795038	8.862	0.4043	1.941	167.848	173.0103	967.9	2.896282
diff --git a/tests/sipnet/test_modeling/termination.param b/tests/sipnet/test_modeling/termination.param
new file mode 100644
index 00000000..fe8fdd45
--- /dev/null
+++ b/tests/sipnet/test_modeling/termination.param
@@ -0,0 +1,78 @@
+plantWoodInit 70.0
+laiInit 2.3481421625603751
+leafCSpWt 29.409172953079818
+litterInit 525.0
+soilInit 5550.0
+soilWFracInit 1.9282707581831169
+snowInit 0
+plantStorageNInit 0.001
+fineRootFrac 0.25
+coarseRootFrac 0.055
+aMax 603.51460591037255
+aMaxFrac 0.75
+psnTMin 2
+psnTOpt 17.796068268198724
+dVpdSlope 0.015793854979347089
+dVpdExp 2
+halfSatPar 15.300351381651099
+attenuation 0.57999999999999996
+baseVegResp 0.0060000000000000001
+baseFolRespFrac 0.016122588146105607
+baseSoilResp 0.030052329942810149
+baseFineRootResp 0.2695450954857716
+baseCoarseRootResp 0.0060000000000000001
+vegRespQ10 1.4019202256757124
+fineRootQ10 2.6000000000000001
+coarseRootQ10 3.20392084014602
+soilRespQ10 1.6759340808079286
+growthRespFrac 0.20000000000000001
+frozenSoilFolREff 0
+frozenSoilThreshold 0
+soilRespMoistEffect 1
+leafOnDay 0
+gddLeafOn 500
+soilTempLeafOn 12
+leafOnReallocFrac 0.10337817506913401
+leafOffDay 0
+leafGrowth 126
+fracLeafFall 1
+woodTurnoverRate 0.014
+leafTurnoverRate 1.09744012449142
+fineRootTurnoverRate 0.2682821481714458
+coarseRootTurnoverRate 0.056000000000000001
+litterBreakdownRate 0.56007424802511896
+fracLitterRespired 0.76768911631115766
+fineRootAllocation 0.10609482754666355
+woodAllocation 0.33470604289674222
+leafAllocation 0.53456185029321135
+waterRemoveFrac 0.087999999999999995
+frozenSoilEff 1
+wueConst 10.9
+soilWHC 12
+immedEvapFrac 0.10000000000000001
+leafPoolDepth 0.10000000000000001
+fastFlowFrac 0
+snowMelt 0.14999999999999999
+rdConst 1149.9961760992301
+rSoilConst1 8.1999999999999993
+rSoilConst2 4.2999999999999998
+cFracLeaf 0.46999999999999997
+mineralNInit 0
+soilOrgNInit 528.09892365915437
+litterOrgNInit 18.860086734898815
+nVolatilizationFrac 6.8860489049188205e-05
+nLeachingFrac 0.248596695361815
+leafNResorptionFrac 0.53599062060008895
+leafCN 15.0
+woodCN 77.0
+fineRootCN 40.0
+kCN 190.0
+nFixationFracMax 0
+halfNFixationMax 0.99630850231037704
+fAnoxia 0.29999999999999999
+anaerobicDecompRate 0.26666853510696198
+anaerobicTransExp 9.9741824548147182
+soilMethaneRate 1.0029297483187848e-05
+litterMethaneRate 4.9804017375032515e-06
+waterDrainFrac 0.10000000000000001
+soilCSaturation 1.0
diff --git a/tests/sipnet/test_modeling/testCompleteHarvest.c b/tests/sipnet/test_modeling/testCompleteHarvest.c
new file mode 100644
index 00000000..d7c3db82
--- /dev/null
+++ b/tests/sipnet/test_modeling/testCompleteHarvest.c
@@ -0,0 +1,560 @@
+// Native three-step fixture from sipnet-termination-repro.tar.gz, derived from
+// wet_d0.1_m1, point 332581. The first step primes mean NPP; positive PAR is
+// intentional even at the midnight harvest timestamp. This is a constructed
+// model regression fixture, not a field observation or restart replay.
+#include "utils/tUtils.h"
+#include "sipnet/events.c"
+#include "sipnet/sipnet.c"
+#include <sys/wait.h>
+
+static int failures;
+static ModelParams *testParams;
+
+static void near(double actual, double expected, const char *name) {
+  if (!isfinite(actual) || fabs(actual - expected) > 1e-8) {
+    logTest("%s: %.15g, expected %.15g\n", name, actual, expected);
+    failures++;
+  }
+}
+
+static void balanced(void) {
+  near(balanceTracker.deltaC, 0, "C balance");
+  near(balanceTracker.deltaN, 0, "N balance");
+  near(balanceTracker.clampedC, 0, "C clamping");
+  near(balanceTracker.clampedN, 0, "N clamping");
+}
+
+static void processEvents(void) {
+  EventNode *event = getCurrentEvent();
+  processEventsForCarbon(event);
+  processEventsForNitrogen(event);
+}
+
+static void start(int mode, const char *harvest) {
+  initContext();
+  ctx.events = 1;
+  ctx.litterPool = mode > 0;
+  ctx.nitrogenCycle = mode == 2;
+  ctx.anaerobic = 1;
+  ctx.waterHResp = 1;
+  ctx.growthResp = 0;
+  ctx.leafWater = 0;
+  ctx.flooding = 1;
+  ctx.gdd = 0;
+  ctx.soilPhenol = 0;
+  envi = (Envi){0};
+  fluxes = (Fluxes){0};
+  eventTrackers = (EventTrackers){0};
+  initModel(&testParams, "termination.param", "termination.clim");
+  setupModel();
+  gEvents = harvest ? createEventNode(2017, 20, HARVEST, harvest) : NULL;
+  setupEvents();
+  openEventOutFile("events.out", 1);
+  updateState();  // Prime mean NPP with the native first timestep.
+  balanced();
+  climate = climate->nextClim;
+}
+
+static void finish(void) {
+  cleanupModel();
+  deleteModelParams(testParams);
+  gEvents = gEvent = NULL;
+}
+
+static void clearPlant(void) {
+  near(envi.plantLeafC, 0, "leaf");
+  near(envi.plantWoodC, 0, "wood");
+  near(envi.fineRootC, 0, "fine roots");
+  near(envi.coarseRootC, 0, "coarse roots");
+  near(envi.plantCAccountingDelta, 0, "accounting C");
+  near(envi.plantStorageN, 0, "storage N");
+  near(plantSurvivalTracker.isAlive, 0, "dead");
+  near(getMeanTrackerMean(meanNPP), 0, "mean NPP");
+}
+
+static void fullCase(int mode, int accounting, int dark, int limited,
+                     const char *harvest, double aboveExport,
+                     double belowExport) {
+  // A no-harvest run supplies the independently calculated end-of-step pools.
+  logTest(
+      "*** Running full case with mode: %d accounting: %d dark: %d limited: %d"
+      " harvest: %s aboveExport: %.2f belowExport: %.2f\n",
+      mode, accounting, dark, limited, harvest ? harvest : "none", aboveExport,
+      belowExport);
+  Envi end = {0};
+  for (int h = 0; h < 2; h++) {
+    start(mode, h ? harvest : NULL);
+    envi.plantCAccountingDelta = accounting;
+    if (dark)
+      climate->par = 0;
+    if (limited && mode == 2) {
+      envi.minN = 0;
+      envi.plantStorageN = 0;
+    }
+    updateState();
+    balanced();
+    if (!h) {
+      // First time through, capture envi state
+      end = envi;
+    } else {
+      const double aboveC =
+          end.plantLeafC + end.plantWoodC + end.plantCAccountingDelta;
+      const double belowC = end.fineRootC + end.coarseRootC;
+      const double aboveN = mode == 2 ? end.plantLeafC / params.leafCN +
+                                            end.plantWoodC / params.woodCN
+                                      : 0;
+      const double belowN = mode == 2 ? end.fineRootC / params.fineRootCN +
+                                            end.coarseRootC / params.woodCN
+                                      : 0;
+      near(envi.soilC,
+           end.soilC + belowC * (1 - belowExport) +
+               (mode == 0 ? aboveC * (1 - aboveExport) : 0),
+           "soil routing");
+      near(envi.litterC, end.litterC + (mode ? aboveC * (1 - aboveExport) : 0),
+           "litter routing");
+      near(fluxes.eventOutputC * climate->length,
+           aboveC * aboveExport + belowC * belowExport, "C export");
+      near(fluxes.eventOutputN * climate->length,
+           aboveN * aboveExport + belowN * belowExport, "N export");
+
+      if (mode == 2) {
+        near(envi.litterN,
+             end.litterN + aboveN * (1 - aboveExport) + end.plantStorageN,
+             "litter N including storage");
+        near(envi.soilOrgN, end.soilOrgN + belowN * (1 - belowExport),
+             "soil N");
+      }
+      clearPlant();
+      climate = climate->nextClim;
+      updateState();
+      balanced();
+      clearPlant();
+      near(fluxes.photosynthesis, 0, "no fallow photosynthesis");
+      // A later planting can restore living tissue.
+      EventNode *plant = createEventNode(2017, 21, PLANTING, "2 3 4 5");
+      gEvents->nextEvent = plant;
+      gEvent = plant;
+      climate->day = 21;
+      updateState();
+      balanced();
+      near(plantSurvivalTracker.isAlive, 1, "alive after planting");
+    }
+    finish();
+  }
+}
+
+static void partialCase(int accounting, const char *harvest, double fraction) {
+  logTest(
+      "*** Running partial case with accounting: %d harvest: %s fraction: %f\n",
+      accounting, harvest, fraction);
+  start(2, harvest);
+  envi.plantCAccountingDelta = accounting;
+  Envi before = envi;
+  resetFluxes();
+  processEvents();
+  updateBalanceTrackerPreUpdate();
+  updatePoolsForEvents();
+  near(envi.plantWoodC, before.plantWoodC * (1 - fraction), "partial wood");
+  near(envi.plantCAccountingDelta, accounting * (1 - fraction),
+       "partial accounting");
+  updateBalanceTrackerPostUpdate();
+  near(balanceTracker.postTotalC + fluxes.eventOutputC * climate->length,
+       balanceTracker.preTotalC, "partial C conservation");
+  near(balanceTracker.postTotalN + fluxes.eventOutputN * climate->length,
+       balanceTracker.preTotalN, "partial N conservation");
+  finish();
+}
+
+static void partialTimestepCase(int accounting) {
+  logTest("*** Running partial timestep case\n");
+  start(2, ".1 .2 .3 .4");
+  envi.plantCAccountingDelta = accounting;
+  updateState();
+  balanced();
+  near(plantSurvivalTracker.isAlive, 1, "partial harvest remains alive");
+  climate = climate->nextClim;
+  updateState();
+  balanced();
+  near(plantSurvivalTracker.isAlive, 1, "partial harvest can continue growth");
+  finish();
+}
+
+static void invalidCase(int which) {
+  pid_t child = fork();
+  if (child < 0) {
+    failures++;
+    return;
+  }
+  if (child == 0) {
+    start(2, "0 0 1 1");
+    if (which < 4) {  // 0, 1, 2, 3
+      const char *bad[] = {"0 0 -1 1", "-0.1 0 1 1", "0 0 1.1 1", "1 0.1 0 1"};
+      // createEventNode should reject these
+      createEventNode(2017, 20, HARVEST, bad[which]);
+    } else {  // 4, 5, 6
+      EventNode *extra = createEventNode(2017, 20, HARVEST, "0 0 .2 .2");
+      if (which == 4) {
+        gEvents->nextEvent = extra;
+        // Harvest greater than 100%
+        processEvents();
+      } else if (which == 5) {
+        extra->nextEvent = gEvents;
+        gEvents = extra;
+        setupEvents();
+        // Harvest greater than 100%
+        processEvents();
+      } else {  // 6
+        EventNode *leafOff1 = createEventNode(2017, 20, LEAFOFF, "");
+        EventNode *leafOff2 = createEventNode(2017, 20, LEAFOFF, "");
+        leafOff1->nextEvent = leafOff2;
+        gEvents = leafOff1;
+        setupEvents();
+        processEvents();
+      }
+    }
+    _exit(99);
+  }
+  int status;
+  waitpid(child, &status, 0);
+  int expected = EXIT_CODE_BAD_PARAMETER_VALUE;
+  if (!WIFEXITED(status) || WEXITSTATUS(status) != expected) {
+    logTest("Invalid case %d did not exit with code %d (was %d)\n", which,
+            expected, WIFEXITED(status) ? WEXITSTATUS(status) : -1);
+    failures++;
+  }
+}
+
+static void coincidentFertilizerCase(void) {
+  logTest("*** Running coincident fertilizer case\n");
+  Envi expected = {0};
+  for (int h = 0; h < 2; h++) {
+    start(2, h ? "0 0 1 1" : NULL);
+    EventNode *fert = createEventNode(2017, 20, FERTILIZATION, "1 2 3");
+    fert->nextEvent = gEvents;
+    gEvents = fert;
+    setupEvents();
+    updateState();
+    balanced();
+    if (!h)
+      expected = envi;
+    else {
+      near(envi.soilC,
+           expected.soilC + expected.fineRootC + expected.coarseRootC,
+           "fertilizer plus harvest soil C");
+      near(envi.litterC,
+           expected.litterC + expected.plantWoodC + expected.plantLeafC +
+               expected.plantCAccountingDelta,
+           "fertilizer plus harvest litter C");
+      near(envi.minN, expected.minN, "fertilizer applied once");
+      clearPlant();
+    }
+    finish();
+  }
+}
+
+static void sameEnvi(const Envi *a, const Envi *b) {
+#define E(f) near(a->f, b->f, #f)
+  E(plantWoodC);
+  E(plantLeafC);
+  E(soilC);
+  E(soilWater);
+  E(litterC);
+  E(snow);
+  E(coarseRootC);
+  E(fineRootC);
+  E(minN);
+  E(soilOrgN);
+  E(litterN);
+  E(plantStorageN);
+  E(plantCAccountingDelta);
+#undef E
+}
+static void coincidentEventCase(int type, const char *arguments) {
+  logTest("*** Running coincident event case with type: %s params: %s\n",
+          eventTypeToString(type), *arguments ? arguments : "<none>");
+  Envi end;
+  Fluxes ordinary;
+  for (int order = 0; order < 3; order++) {
+    start(2, order ? ".25 .75 .75 .25" : NULL);
+    params.leafGrowth = 1;
+    params.fracLeafFall = .25;
+    EventNode *extra = createEventNode(2017, 20, type, arguments);
+    if (order == 2)
+      gEvents->nextEvent = extra;
+    else {
+      extra->nextEvent = gEvents;
+      gEvents = extra;
+    }
+    setupEvents();
+    updateState();
+    balanced();
+    if (!order) {
+      end = envi;
+      ordinary = fluxes;
+      if ((type == LEAFON && fluxes.eventLeafOnCreation <= 0) ||
+          (type == LEAFOFF && fluxes.eventLeafOffLitterC <= 0) ||
+          (type == IRRIGATION && fluxes.eventSoilWater <= 0))
+        failures++;
+    } else {
+      double dt = climate->length;
+      double ac = end.plantWoodC + end.plantCAccountingDelta + end.plantLeafC;
+      double bc = end.fineRootC + end.coarseRootC;
+      double an =
+          end.plantWoodC / params.woodCN + end.plantLeafC / params.leafCN;
+      double bn =
+          end.fineRootC / params.fineRootCN + end.coarseRootC / params.woodCN;
+      Envi expected = end;
+      expected.soilC += .25 * bc;
+      expected.litterC += .75 * ac;
+      expected.soilOrgN += .25 * bn;
+      expected.litterN += .75 * an + end.plantStorageN;
+      expected.plantWoodC = expected.plantLeafC = expected.fineRootC = 0;
+      expected.coarseRootC = expected.plantCAccountingDelta =
+          expected.plantStorageN = 0;
+      sameEnvi(&envi, &expected);
+      clearPlant();
+      near(fluxes.eventInputC, ordinary.eventInputC, "C input once");
+      near(fluxes.eventInputN, ordinary.eventInputN, "N input once");
+      near(fluxes.eventOutputC,
+           ordinary.eventOutputC + (.25 * ac + .75 * bc) / dt, "C export");
+      near(fluxes.eventOutputN,
+           ordinary.eventOutputN + (.25 * an + .75 * bn) / dt, "N export");
+#define F(f) near(fluxes.f, ordinary.f, #f)
+      F(eventSoilWater);
+      F(eventEvap);
+      F(eventLeafOnCreation);
+      F(eventLeafOnCreationFromWood);
+      F(eventLeafOffLitterC);
+      F(eventLeafOffNResorption);
+#undef F
+      climate = climate->nextClim;
+      updateState();
+      balanced();
+      clearPlant();
+      near(fluxes.photosynthesis, 0, "no fallow photosynthesis");
+    }
+    finish();
+  }
+}
+
+static void restartCase(void) {
+  logTest("*** Running restart case\n");
+  Envi expected;
+  for (int resumed = 0; resumed < 2; resumed++) {
+    start(2, "0 0 1 1");
+    if (resumed) {
+      restartNoteProcessedClimateStep(firstClimate);
+      restartWriteCheckpoint("termination.restart", meanNPP);
+      envi = (Envi){0};
+      resetMeanTracker(meanNPP, 0);
+      restartLoadCheckpoint("termination.restart", meanNPP);
+    }
+    updateState();
+    balanced();
+    clearPlant();
+    if (resumed) {
+      restartNoteProcessedClimateStep(climate);
+      restartWriteCheckpoint("termination.restart", meanNPP);
+    }
+    climate = climate->nextClim;
+    if (resumed) {
+      envi = (Envi){0};
+      restartLoadCheckpoint("termination.restart", meanNPP);
+    }
+    updateState();
+    balanced();
+    clearPlant();
+    if (!resumed)
+      expected = envi;
+    else {
+      sameEnvi(&envi, &expected);
+      near(envi.soilC, expected.soilC, "restart soil C");
+      near(envi.litterC, expected.litterC, "restart litter C");
+      near(envi.minN, expected.minN, "restart mineral N");
+      near(envi.soilOrgN, expected.soilOrgN, "restart soil N");
+      near(envi.litterN, expected.litterN, "restart litter N");
+    }
+    finish();
+  }
+}
+
+static void leafBudgetCases(void) {
+  int caseNum = 0;
+  char harvests[4][16] = {"none", "0 0 1 1", "0.5 0.5 0.5 0.5", "1 1 0 0"};
+  for (int kind = 0; kind < 4; kind++) {
+    // kind:          0    1    2    3
+    // fracLeafFall   1  .25  .40  .75
+    // leafOff        a    a    b    c
+    // a: inserted at start
+    // b: two events inserted at start
+    // c: no leaf-off events
+    for (int dark = 0; dark < 2; dark++) {
+      // dark : set climate->par=0 when dark
+      for (int sign = -1; sign <= 1; sign++) {
+        // sign: init npp tracker with 20*sign (-20, 0, 20)
+        for (int account = -1; account <= 1; account++) {
+          // account: starting value for plantCAccountingDelta
+          for (int limited = 0; limited < 2; limited++) {
+            // limited: starting minN = 0 if limited (1000 else)
+            for (int resorb = 0; resorb < 3; resorb++) {
+              /// resorb: leafNResorptionFrac = 0, .5, 1  (resorb/2)
+              for (int harvest = 0; harvest < 4; harvest++) {
+                char *harvestStr = harvests[harvest];
+                logTest(
+                    "*** Running leaf budget case [%d] with kind: %d dark: %d "
+                    "sign %d account %d limited %d resorb %.1f harvest %s\n",
+                    caseNum++, kind, dark, sign, account, limited, resorb / 2.0,
+                    harvestStr);
+                start(2, harvest ? harvestStr : NULL);
+                envi.plantLeafC = 1;
+                envi.plantWoodC = envi.fineRootC = envi.coarseRootC = 100;
+                envi.plantCAccountingDelta = account;
+                envi.plantStorageN = 0;
+                envi.minN = limited ? 0 : 1000;
+                params.leafAllocation = params.woodAllocation = .25;
+                params.fineRootAllocation = params.coarseRootAllocation = .25;
+                params.leafTurnoverRate = 2;
+                params.woodTurnoverRate = params.fineRootTurnoverRate = 0;
+                params.coarseRootTurnoverRate = params.litterBreakdownRate = 0;
+                params.baseSoilResp = params.soilMethaneRate =
+                    params.litterMethaneRate = 0;
+                params.nVolatilizationFrac = params.nLeachingFrac = 0;
+                params.nFixationFracMax = 0;
+                params.leafNResorptionFrac = .5 * resorb;
+                params.fracLeafFall = kind == 0   ? 1
+                                      : kind == 1 ? .25
+                                      : kind == 2 ? .40
+                                                  : .75;
+                if (dark)
+                  climate->par = 0;
+                resetMeanTracker(meanNPP, 20 * sign);
+                if (kind != 3) {
+                  EventNode *off = createEventNode(2017, 20, LEAFOFF, "");
+                  off->nextEvent = gEvents;
+                  gEvents = off;
+                  if (kind == 2) {
+                    EventNode *second = createEventNode(2017, 20, LEAFOFF, "");
+                    second->nextEvent = off->nextEvent;
+                    off->nextEvent = second;
+                  }
+                  setupEvents();
+                }
+                double beforeLeafC = envi.plantLeafC;
+                double beforeLitter = envi.litterC;
+                double beforeLitterN = envi.litterN;
+                updateState();
+                balanced();
+                double shed = (fluxes.leafLitter + fluxes.leafOffLitter +
+                               fluxes.eventLeafOffLitterC) *
+                              climate->length;
+                double pool = 1 + fluxes.leafCreation * climate->length;
+                if (pool - shed < -1e-10 || envi.plantLeafC < -1e-10) {
+                  logTest("Leaf budget overdraw (2): shed %.17g pool %.17g "
+                          "|s-p| %.17g leaf %.17g expr1 %d expr2 %d\n",
+                          shed, pool, pool - shed, envi.plantLeafC,
+                          pool - shed < -1e-10, envi.plantLeafC < -1e-10);
+                  failures++;
+                }
+
+                double expLeafOffLitter =
+                    kind == 3 ? 0
+                              : fmin(params.fracLeafFall * (1 + (kind == 2)),
+                                     beforeLeafC + (fluxes.leafCreation -
+                                                    fluxes.leafLitter) *
+                                                       climate->length);
+                near(fluxes.eventLeafOffLitterC * climate->length,
+                     expLeafOffLitter, "event leaves conserved");
+                if (!harvest) {
+                  near(envi.litterC - beforeLitter, shed,
+                       "leaf litter transfer");
+                  near(envi.litterN - beforeLitterN,
+                       shed * (1 - params.leafNResorptionFrac) / params.leafCN,
+                       "leaf litter N transfer");
+                } else {
+                  clearPlant();
+                  climate = climate->nextClim;
+                  updateState();
+                  balanced();
+                  clearPlant();
+                  near(fluxes.photosynthesis, 0,
+                       "no regrowth after complete harvest");
+                }
+                finish();
+              }  // harvest loop
+            }  // resorb loop
+          }  // limited loop
+        }  // account loop
+      }  // sign loop
+    }  // dark loop
+  }  // kind loop
+}
+
+// Full export during net loss must not borrow carbon from an empty soil pool.
+static void exportLossCase(void) {
+  start(1, "1 1 0 0");
+  envi.litterC = envi.soilC = envi.litterN = envi.soilOrgN = 0;
+  params.baseSoilResp = params.litterBreakdownRate = 0;
+  params.soilMethaneRate = params.litterMethaneRate = 0;
+  climate->par = 0;
+  resetMeanTracker(meanNPP, -20);
+  updateState();
+  balanced();
+  clearPlant();
+  finish();
+}
+
+int main(void) {
+
+  logTest("Starting testCompleteHarvest\n");
+
+  const char *harvest[] = {"0 0 1 1", "1 1 0 0", ".25 .75 .75 .25"};
+  const double above[] = {0, 1, .25}, below[] = {0, 1, .75};
+  for (int mode = 0; mode < 3; mode++)
+    for (int account = -1; account <= 1; account++)
+      for (int dark = 0; dark < 2; dark++)
+        for (int limited = 0; limited <= (mode == 2); limited++)
+          for (int route = 0; route < 3; route++)
+            fullCase(mode, account, dark, limited, harvest[route], above[route],
+                     below[route]);
+
+  logTest("\n");
+  for (int account = -1; account <= 1; account++) {
+    partialTimestepCase(account);
+    partialCase(account, ".1 .2 .3 .4", .4);
+    partialCase(account, "0 0 .999999999 .999999999", .999999999);
+  }
+
+  logTest("\n");
+  coincidentFertilizerCase();
+  coincidentEventCase(LEAFON, "");
+  coincidentEventCase(LEAFOFF, "");
+  coincidentEventCase(IRRIGATION, "2 0");
+  restartCase();
+
+  logTest("\n");
+
+  leafBudgetCases();
+  exportLossCase();
+
+  // Leave these last, as the forking messes up the output otherwise
+  logTest("\n");
+  logTest("*** Running invalid case checks; seven errors expected\n");
+  for (int which = 0; which < 7; which++) {
+    invalidCase(which);
+  }
+
+  logTest("\n");
+  logTest("Complete harvest total failures: %d\n", failures);
+  logTest("\n");
+
+  int status = failures > 0;
+
+  if (status) {
+    logTest("FAILED testCompleteHarvest with status %d\n", status);
+    exit(status);
+  }
+
+  logTest("PASSED testCompleteHarvest\n");
+
+  return failures != 0;
+}
diff --git a/tests/sipnet/test_modeling/testNitrogenCycle.c b/tests/sipnet/test_modeling/testNitrogenCycle.c
index 95e419c0..1e061e51 100644
--- a/tests/sipnet/test_modeling/testNitrogenCycle.c
+++ b/tests/sipnet/test_modeling/testNitrogenCycle.c
@@ -2,6 +2,7 @@
 #include "sipnet/events.c"
 #include "sipnet/nitrogen.c"
 #include "sipnet/limitations.c"
+#include "utils/helpers.c"
 
 /////
 // Setup and general test state management
@@ -134,15 +135,17 @@ int testFertilization(void) {
   // init minN 2, nVol 0.1
   initNVolatilizationState(initN, nVolFrac);
 
-  // fert event: 15 5 10
+  // fert event: 15 orgN 5 orgC 10 minN
   double fertMinN = 10;
   initEvents("events_fert.in", "events.out", 0);
   setupEvents();
 
+  EventNode *event = getCurrentEvent();
+  processEventsForCarbon(event);
+  processEventsForNitrogen(event);
   calcNVolatilizationFlux();
-  processEvents();
-  updateNitrogenPools();
   updatePoolsForEvents();
+  updateNitrogenPools();
 
   // Want to test:
   // envi: minN
@@ -708,6 +711,8 @@ int testLeafTurnoverNResorption(void) {
   envi.minN = 1.0;
   envi.plantStorageN = 0.0;
   fluxes.leafOffNResorption = 2.0;
+  params.leafNResorptionFrac = 0.5;
+
   updateNitrogenPools();
 
   double expStorageN = 2.0 * climate->length;  // 0.25
@@ -727,9 +732,14 @@ int testLeafTurnoverNResorption(void) {
   double availableN = minN2;  // plus unclaimed plantStorageN, which is 0 here
   double maxUptake = demandFlux * climate->length;
   double reduction = availableN / maxUptake;
+  // enough so we don't limit leaf litter in checkNegativeCreation
+  envi.plantLeafC = 8.0;
 
   initNLimitationState(minN2, 0);
   fluxes.leafOffNResorption = resorpFlux;
+  // to make sure leafOffNResorption is not reduced by leafOffLitter
+  fluxes.leafOffLitter =
+      resorpFlux * params.leafCN / params.leafNResorptionFrac;
 
   doNFixUpLimitCalcs();
   updateNitrogenPools();
@@ -738,7 +748,7 @@ int testLeafTurnoverNResorption(void) {
                                  "[turnover resorption] leafCreation");
   status |= checkNLimitationFlux(fluxes.woodCreation, 500 * reduction,
                                  "[turnover resorption] woodCreation");
-  status |= checkMinAndStorageN("turnover resorption", 0.0,
+  status |= checkMinAndStorageN("turnover resorption final", 0.0,
                                 resorpFlux * climate->length);
 
   return status;
diff --git a/tests/sipnet/test_modeling/testPlantMortality.c b/tests/sipnet/test_modeling/testPlantMortality.c
index 9e5d3b36..d7bace6d 100644
--- a/tests/sipnet/test_modeling/testPlantMortality.c
+++ b/tests/sipnet/test_modeling/testPlantMortality.c
@@ -44,8 +44,8 @@ void resetEnv(void) {
   envi.soilOrgN = 0.0;
   envi.litterN = 0.0;
   envi.plantStorageN = 0.0;
-  eventTrackers.harvestFracRemoved = 0.0;
-  eventTrackers.harvestFracTransferred = 0.0;
+  eventTrackers.harvestTrackers.totalFracRemoved = 0.0;
+  eventTrackers.harvestTrackers.totalFracTransferred = 0.0;
 }
 
 int checkPool(double calc, double exp, const char *label) {
diff --git a/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c b/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c
index 793712b2..d5d918a3 100644
--- a/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c
+++ b/tests/sipnet/test_sipnet_infrastructure/testDebugLogFiles.c
@@ -3,6 +3,7 @@
 #include <string.h>
 
 #include "common/logging.h"
+#include "sipnet/events.h"
 #include "sipnet/state.h"
 #include "utils/tUtils.h"
 
@@ -15,14 +16,12 @@
 
 #define EXPECTED_ENVI_FIELDS ((int)(sizeof(Envi) / sizeof(double)))
 #define EXPECTED_FLUX_FIELDS ((int)(sizeof(Fluxes) / sizeof(double)))
-// These are not all doubles
-// Let's do a little more work for Trackers
-#define NUM_TRACKER_INTS 1
-#define TRACKERS_SIZE (int)(sizeof(Trackers) - NUM_TRACKER_INTS * sizeof(int))
-#define EXPECTED_TRACKER_FIELDS                                                \
-  ((int)(TRACKERS_SIZE / sizeof(double)) + NUM_TRACKER_INTS)
+// The Trackers struct is not all doubles, but the int(s) get padded to 8 bytes
+#define EXPECTED_TRACKER_FIELDS ((int)(sizeof(Trackers) / sizeof(double)))
 #define EXPECTED_PHENOLOGY_FIELDS 3
 #define EXPECTED_SURVIVAL_FIELDS 1
+#define EXPECTED_EVENT_TRACKER_FIELDS                                          \
+  (int)(sizeof(EventTrackers) / sizeof(double))
 
 static int countTokens(const char *line) {
   int count = 0;
@@ -114,7 +113,7 @@ int run(void) {
   status |=
       checkHeader(TRACKERS_FILE,
                   3 + EXPECTED_TRACKER_FIELDS + EXPECTED_PHENOLOGY_FIELDS +
-                      EXPECTED_SURVIVAL_FIELDS,
+                      EXPECTED_SURVIVAL_FIELDS + EXPECTED_EVENT_TRACKER_FIELDS,
                   "t.gpp", "pt.lastYear");
 
   int mainLines = countLines(SIPNET_OUT_FILE);
diff --git a/tests/smoke/russell_2/events.out b/tests/smoke/russell_2/events.out
index 45b9422a..f387a652 100644
--- a/tests/smoke/russell_2/events.out
+++ b/tests/smoke/russell_2/events.out
@@ -1,11 +1,11 @@
 year  day  type     param_name=delta[,param_name=delta,...]
 2016   47  leafon   leafOnCreation=114.61,leafOnCreationFromWood=86.11
-2016   90  fert     eventLitterC=5.00,eventSoilC=0.00,eventMinN=10.00,eventLitterN=15.00,eventInputC=5.00,eventInputN=25.00
+2016   90  fert     eventLitterC=5.00,eventSoilC=0.00,eventInputC=5.00,eventMinN=10.00,eventLitterN=15.00,eventInputN=25.00
 2016   92  irrig    eventSoilWater=2.80,eventEvap=0.00
 2016   96  irrig    eventSoilWater=2.80,eventEvap=0.00
 2016  100  irrig    eventSoilWater=2.80,eventEvap=0.00
 2016  104  irrig    eventSoilWater=2.80,eventEvap=0.00
-2016  104  fert     eventLitterC=5.00,eventSoilC=0.00,eventMinN=10.00,eventLitterN=15.00,eventInputC=5.00,eventInputN=25.00
+2016  104  fert     eventLitterC=5.00,eventSoilC=0.00,eventInputC=5.00,eventMinN=10.00,eventLitterN=15.00,eventInputN=25.00
 2016  108  irrig    eventSoilWater=2.80,eventEvap=0.00
 2016  112  irrig    eventSoilWater=2.80,eventEvap=0.00
 2016  116  irrig    eventSoilWater=2.80,eventEvap=0.00
diff --git a/tests/utils/helpers.c b/tests/utils/helpers.c
index ff2e5b91..01fc1959 100644
--- a/tests/utils/helpers.c
+++ b/tests/utils/helpers.c
@@ -21,6 +21,9 @@ void prepTypesTest(void) {
 
 void procEvents(void) {
   resetFluxes();
-  processEvents();
+  EventNode *event = getCurrentEvent();
+  processEventsForCarbon(event);
+  processEventsForNitrogen(event);
+  writeEventsOut();
   updatePoolsForEvents();
 }