diff --git a/include/topmodel.h b/include/topmodel.h
index 6fe0306..cb99686 100755
--- a/include/topmodel.h
+++ b/include/topmodel.h
@@ -68,7 +68,7 @@ extern void topmod(FILE *output_fptr, int nstep, int num_topodex_values,
                 const double *time_delay_histogram,
                 double *sbar,int num_delay, int current_time_step, int stand_alone,
                 double *sump, double *sumae, double *sumq, double *sumrz, double *sumuz,
-                double *quz, double *qb, double *qof, double *p, double *ep);
+                double *quz, double *qb, double *qof, double *p, double *ep, double *ponded_depth);
 
 extern int tread(FILE *subcat_fptr,FILE *output_fptr,char *subcat,
                 int *num_topodex_values,int *num_channels,double *area,
@@ -184,6 +184,7 @@ struct TopModel_Struct{
   double qof; /* flow from saturated area and infiltration excess flow*/
   double p;   /* adjusted rain*/
   double ep;  /* adjusted potential evaporation*/
+  double ponded_depth; /* queued delayed flow from hydrograph ordinates */
 
   /************** Framework vars **************/ 
   int stand_alone;
diff --git a/src/bmi_serialization.cpp b/src/bmi_serialization.cpp
index 4daa17d..9647f60 100644
--- a/src/bmi_serialization.cpp
+++ b/src/bmi_serialization.cpp
@@ -54,6 +54,7 @@ void TopmodelSerializer::serialize(Archive& ar, const unsigned int version) {
     ar & model->qof; //
     ar & model->p; // reassigned each update; used in calc after assignment
     ar & model->ep; // reassigned each update; used in calc after assignment
+    ar & model->ponded_depth; // reassigned each update
     ar & model->sbar; // used then reassigned
 
     // array data that updates in update; counts set in config
diff --git a/src/bmi_topmodel.c b/src/bmi_topmodel.c
index ef39709..ec8a9e1 100755
--- a/src/bmi_topmodel.c
+++ b/src/bmi_topmodel.c
@@ -6,7 +6,7 @@
 
 /* BMI Adaption: Max i/o file name length changed from 30 to 256 */
 #define MAX_FILENAME_LENGTH   256
-#define OUTPUT_VAR_NAME_COUNT 14
+#define OUTPUT_VAR_NAME_COUNT 15
 #define INPUT_VAR_NAME_COUNT  2
 #define PARAM_VAR_NAME_COUNT  8
 
@@ -27,7 +27,8 @@ static const char *output_var_names[OUTPUT_VAR_NAME_COUNT] = {
     "land_surface_water__domain_time_integral_of_runoff_volume_flux", // sumq
     "soil_water__domain_root-zone_volume_deficit", // sumrz
     "soil_water__domain_unsaturated-zone_volume", // sumuz
-    "land_surface_water__water_balance_volume" // bal
+    "land_surface_water__water_balance_volume", // bal
+    "nwm_ponded_depth" // sum of Q[1..num_time_delay_histo_ords]
 };
 
 static const char *output_var_types[OUTPUT_VAR_NAME_COUNT] = {
@@ -44,11 +45,12 @@ static const char *output_var_types[OUTPUT_VAR_NAME_COUNT] = {
     "double",
     "double",
     "double",
+    "double",
     "double"
 };
 
 static const int output_var_item_count[OUTPUT_VAR_NAME_COUNT] =
-    {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+    {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
 
 static const char *output_var_units[OUTPUT_VAR_NAME_COUNT] = {
     "m h-1",
@@ -64,11 +66,12 @@ static const char *output_var_units[OUTPUT_VAR_NAME_COUNT] = {
     "m",
     "m",
     "m",
+    "m",
     "m"
 };
 
 static const int output_var_grids[OUTPUT_VAR_NAME_COUNT] =
-    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
 static const char *output_var_locations[OUTPUT_VAR_NAME_COUNT] = {
     "node",
@@ -84,6 +87,7 @@ static const char *output_var_locations[OUTPUT_VAR_NAME_COUNT] = {
     "node",
     "node",
     "node",
+    "node",
     "node"
 };
 
@@ -501,7 +505,8 @@ static int Update(Bmi *self) {
         &topmodel->qb,
         &topmodel->qof,
         &topmodel->p,
-        &topmodel->ep
+        &topmodel->ep,
+	&topmodel->ponded_depth
     );
 
     return BMI_SUCCESS;
@@ -911,6 +916,14 @@ static int Get_value_ptr(Bmi *self, const char *name, void **dest) {
         *dest    = (void *)&topmodel->bal;
         return BMI_SUCCESS;
     }
+
+    // ponded depth
+    if (strcmp(name, "nwm_ponded_depth") == 0) {
+        topmodel_model *topmodel;
+        topmodel = (topmodel_model *)self->data;
+        *dest    = (void *)&topmodel->ponded_depth;
+        return BMI_SUCCESS;
+    }
     // szm (parameter)
     if (strcmp(name, "szm") == 0) {
         topmodel_model *topmodel;
diff --git a/src/topmodel.c b/src/topmodel.c
index 4eeea89..9cc3cc0 100755
--- a/src/topmodel.c
+++ b/src/topmodel.c
@@ -147,7 +147,8 @@ extern void topmod(
     double *qb,
     double *qof,
     double *p,
-    double *ep
+    double *ep,
+    double *ponded_depth
 ) {
     /*****************************************************************
 
@@ -357,6 +358,13 @@ extern void topmod(
         Q[in] += (*Qout) * time_delay_histogram[ir];
     }
 
+    /* Ponded depth: sum delayed-flow components currently stored in the
+     * hydrograph ordinates portion of the routing array. */
+    *ponded_depth = 0.0;
+    for (ir = 1; ir <= num_time_delay_histo_ords; ir++) {
+        *ponded_depth += Q[ir];
+    }
+
     // Add current time flow to mass balance variable
     *sumq += Q[it];
     /* BMI Adaption: replace nstep with current_time_step */