diff --git a/R/RcppExports.R b/R/RcppExports.R
index 93520faaa..3d6a713c9 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -3293,6 +3293,26 @@ FF16_Environment__ctor <- function(canopy_rescale_usually, soil_number_of_depths
     .Call('_plant_FF16_Environment__ctor', PACKAGE = 'plant', canopy_rescale_usually, soil_number_of_depths)
 }
 
+FF16_Environment__set_extrinsic_driver <- function(obj_, driver_name, x, y) {
+    invisible(.Call('_plant_FF16_Environment__set_extrinsic_driver', PACKAGE = 'plant', obj_, driver_name, x, y))
+}
+
+FF16_Environment__extrinsic_driver_extrapolation <- function(obj_, driver_name, extrapolate) {
+    invisible(.Call('_plant_FF16_Environment__extrinsic_driver_extrapolation', PACKAGE = 'plant', obj_, driver_name, extrapolate))
+}
+
+FF16_Environment__extrinsic_driver_evaluate <- function(obj_, driver_name, u) {
+    .Call('_plant_FF16_Environment__extrinsic_driver_evaluate', PACKAGE = 'plant', obj_, driver_name, u)
+}
+
+FF16_Environment__extrinsic_driver_evaluate_range <- function(obj_, driver_name, u) {
+    .Call('_plant_FF16_Environment__extrinsic_driver_evaluate_range', PACKAGE = 'plant', obj_, driver_name, u)
+}
+
+FF16_Environment__get_extrinsic_driver_names <- function(obj_) {
+    .Call('_plant_FF16_Environment__get_extrinsic_driver_names', PACKAGE = 'plant', obj_)
+}
+
 FF16_Environment__canopy_openness <- function(obj_, height) {
     .Call('_plant_FF16_Environment__canopy_openness', PACKAGE = 'plant', obj_, height)
 }
@@ -3305,10 +3325,6 @@ FF16_Environment__set_fixed_environment <- function(obj_, value, height_max) {
     invisible(.Call('_plant_FF16_Environment__set_fixed_environment', PACKAGE = 'plant', obj_, value, height_max))
 }
 
-FF16_Environment__set_soil_infiltration_rate <- function(obj_, rate) {
-    invisible(.Call('_plant_FF16_Environment__set_soil_infiltration_rate', PACKAGE = 'plant', obj_, rate))
-}
-
 FF16_Environment__set_soil_water_state <- function(obj_, state) {
     invisible(.Call('_plant_FF16_Environment__set_soil_water_state', PACKAGE = 'plant', obj_, state))
 }
diff --git a/R/RcppR6.R b/R/RcppR6.R
index 67d57ccdc..7c5a8a115 100644
--- a/R/RcppR6.R
+++ b/R/RcppR6.R
@@ -1,6 +1,6 @@
 ## Generated by RcppR6: do not edit by hand
 ## Version: 0.2.4
-## Hash: f6b511a10ca3bb03a599e71a6fa83092
+## Hash: 2cb4fee1111abff105bb05242a8b85a0
 
 ##' @importFrom Rcpp evalCpp
 ##' @importFrom R6 R6Class
@@ -4814,6 +4814,21 @@ StochasticPatchRunner <- function(T, E) {
       initialize = function(ptr) {
         self$.ptr <- ptr
       },
+      set_extrinsic_driver = function(driver_name, x, y) {
+        FF16_Environment__set_extrinsic_driver(self, driver_name, x, y)
+      },
+      extrinsic_driver_extrapolation = function(driver_name, extrapolate) {
+        FF16_Environment__extrinsic_driver_extrapolation(self, driver_name, extrapolate)
+      },
+      extrinsic_driver_evaluate = function(driver_name, u) {
+        FF16_Environment__extrinsic_driver_evaluate(self, driver_name, u)
+      },
+      extrinsic_driver_evaluate_range = function(driver_name, u) {
+        FF16_Environment__extrinsic_driver_evaluate_range(self, driver_name, u)
+      },
+      get_extrinsic_driver_names = function() {
+        FF16_Environment__get_extrinsic_driver_names(self)
+      },
       canopy_openness = function(height) {
         FF16_Environment__canopy_openness(self, height)
       },
@@ -4823,9 +4838,6 @@ StochasticPatchRunner <- function(T, E) {
       set_fixed_environment = function(value, height_max) {
         FF16_Environment__set_fixed_environment(self, value, height_max)
       },
-      set_soil_infiltration_rate = function(rate) {
-        FF16_Environment__set_soil_infiltration_rate(self, rate)
-      },
       set_soil_water_state = function(state) {
         FF16_Environment__set_soil_water_state(self, state)
       },
diff --git a/R/ff16w.R b/R/ff16w.R
index fc089ef79..6e85b1fc8 100644
--- a/R/ff16w.R
+++ b/R/ff16w.R
@@ -64,16 +64,15 @@ FF16w_StochasticPatchRunner <- function(p) {
 ## Helper:
 ##' @export
 ##' @rdname FF16_Environment
-##' @param infil_rate rate of water entering the first layer
-##' @param n_layers the number of layers
-##' @param init starting conditions
+##' @param soil_number_of_depths the number of soil layers
+##' @param rainfall constant function value for rainfall driver, y = rainfall
 FF16w_make_environment <- function(canopy_light_tol = 1e-4, 
                                    canopy_light_nbase = 17,
                                    canopy_light_max_depth = 16, 
                                    canopy_rescale_usually = TRUE,
                                    soil_number_of_depths = 1,
                                    soil_initial_state = 0.0,
-                                   soil_infiltration_rate = 0.0) {
+                                   rainfall = 1) {
   
   if(soil_number_of_depths < 1)
     stop("FF16w Environment must have at least one soil layer")
@@ -85,7 +84,7 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4,
                      canopy_light_nbase, 
                      canopy_light_max_depth)
 
-  # there might be a better way to skip this if using defaultss
+  # there might be a better way to skip this if using defaults
   if(sum(soil_initial_state) > 0.0) {
     if(soil_number_of_depths != length(soil_initial_state))
       stop("Not enough starting points for all layers")
@@ -93,7 +92,10 @@ FF16w_make_environment <- function(canopy_light_tol = 1e-4,
     e$set_soil_water_state(soil_initial_state)
   }
     
-  e$set_soil_infiltration_rate(soil_infiltration_rate)
+  x <- seq(0, 10, 1)
+  y <- rep(rainfall, 11)
+  e$set_extrinsic_driver("rainfall", x, y)
+  e$extrinsic_driver_extrapolation("rainfall", TRUE);
   
   return(e)
 }
diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml
index 815517dbe..2c9a50572 100644
--- a/inst/RcppR6_classes.yml
+++ b/inst/RcppR6_classes.yml
@@ -871,6 +871,22 @@ FF16_Environment:
       FF16_Environment object
       @export
   methods:
+    # unforunately need to expose these base Environment functions for every derived class...
+    set_extrinsic_driver:
+      args: [driver_name: "std::string", x: "std::vector<double>", y: "std::vector<double>"]
+      return_type: void
+    extrinsic_driver_extrapolation:
+      args: [driver_name: "std::string", extrapolate: bool]
+      return_type: void
+    extrinsic_driver_evaluate:
+      args: [driver_name: "std::string", u: double]
+      return_type: double
+    extrinsic_driver_evaluate_range:
+      args: [driver_name: "std::string", u: "std::vector<double>"]
+      return_type: "std::vector<double>"
+    get_extrinsic_driver_names:
+      args: []
+      return_type: "std::vector<std::string>"
     canopy_openness:
       args: [height: double]
       return_type: double
@@ -879,9 +895,6 @@ FF16_Environment:
     set_fixed_environment:
       args: [value: double, height_max: double]
       return_type: void
-    set_soil_infiltration_rate:
-      args: [rate: double]
-      return_type: void
     set_soil_water_state:
       args: [state: "std::vector<double>"]
       return_type: void
diff --git a/inst/include/plant/environment.h b/inst/include/plant/environment.h
index 67073917d..c5edea29f 100644
--- a/inst/include/plant/environment.h
+++ b/inst/include/plant/environment.h
@@ -8,6 +8,7 @@
 #include <plant/ode_interface.h>
 #include <plant/internals.h>
 #include <plant/util.h>
+#include <unordered_map>
 #include <Rcpp.h>
 
 using namespace Rcpp;
@@ -69,6 +70,44 @@ class Environment {
   size_t species_arriving_index;
 
   Internals vars;
+
+  /* EXSTRINSIC DRIVERS API*/
+
+  // initialise spline of driver with x, y control points
+  void set_extrinsic_driver(std::string driver_name, std::vector<double> const& x, std::vector<double> const& y) {
+    // if we wanted to be faster we could skip this check (but less safe)
+    if (extrinsic_drivers.count(driver_name) == 0) {
+      util::stop(driver_name + " doesn't exist in the list of extrinsic_drivers.");
+    } else {
+      extrinsic_drivers.at(driver_name).init(x, y);
+      extrinsic_driver_extrapolation(driver_name, false); // default no extrapolation
+    }
+  }
+
+  void extrinsic_driver_extrapolation(std::string driver_name, bool extrapolate) {
+    extrinsic_drivers.at(driver_name).setExtrapolate(extrapolate);
+  }
+
+  // evaluate/query interpolated spline for driver at point u, return s(u), where s is interpolated function
+  double extrinsic_driver_evaluate(std::string driver_name, double u) const {
+    return extrinsic_drivers.at(driver_name).eval(u);
+  }
+
+  // // evaluate/query interpolated spline for driver at vector of points, return vector of values
+  std::vector<double> extrinsic_driver_evaluate_range(std::string driver_name, std::vector<double> u) const {
+    return extrinsic_drivers.at(driver_name).r_eval(u);
+  }
+
+  std::vector<std::string> get_extrinsic_driver_names() {
+    auto ret = std::vector<std::string>();
+    for (auto const& driver : extrinsic_drivers) {
+      ret.push_back(driver.first);
+    }
+    return ret;
+  }
+
+protected:
+  std::unordered_map<std::string, interpolator::Interpolator> extrinsic_drivers;
 };
 
 }
diff --git a/inst/include/plant/interpolator.h b/inst/include/plant/interpolator.h
index 4fca2f521..738e2de16 100644
--- a/inst/include/plant/interpolator.h
+++ b/inst/include/plant/interpolator.h
@@ -16,16 +16,20 @@ class Interpolator {
   void initialise();
 
   void add_point(double xi, double yi);
+  void add_point_sorted(double xi, double yi);
   void clear();
 
   double eval(double u) const;
+  double operator()(double u) const;
   size_t size() const;
 
   double min() const;
   double max() const;
+  void setExtrapolate(bool e);
 
   std::vector<double> get_x() const;
   std::vector<double> get_y() const;
+  
 
   // * R interface
   SEXP r_get_xy() const;
@@ -35,7 +39,8 @@ class Interpolator {
   void check_active() const;
   std::vector<double> x, y;
   tk::spline tk_spline;
-  bool active;
+  bool active = false;
+  bool extrapolate = true;
 };
 
 }
diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h
index 0b2b8cd0f..563e84240 100644
--- a/inst/include/plant/models/ff16_environment.h
+++ b/inst/include/plant/models/ff16_environment.h
@@ -4,6 +4,8 @@
 
 #include <plant/environment.h>
 #include <plant/canopy.h>
+#include <plant/interpolator.h>
+#include <iostream>
 
 using namespace Rcpp;
 
@@ -19,14 +21,16 @@ class FF16_Environment : public Environment {
       : canopy_rescale_usually(canopy_rescale_usually) {
     time = 0.0;
     canopy = Canopy();
-    soil_infiltration_rate = 0.0;
     vars = Internals(soil_number_of_depths);
     set_soil_water_state(std::vector<double>(soil_number_of_depths, 0.0));
+    extrinsic_drivers["rainfall"] = interpolator::Interpolator();
   };
 
+
   // Light interface
   bool canopy_rescale_usually;
 
+  // private?
   Canopy canopy;
 
   // Should this be here or in canopy?
@@ -51,12 +55,9 @@ class FF16_Environment : public Environment {
     canopy.r_init_interpolators(state);
   }
 
-  // Soil interface
-  double soil_infiltration_rate;
-
   virtual void compute_rates() {
     for (size_t i = 0; i < vars.state_size; i++) {
-      vars.set_rate(i, soil_infiltration_rate / (i+1));
+      vars.set_rate(i, extrinsic_driver_evaluate("rainfall", time) / (i+1));
     }
   }
 
@@ -74,10 +75,6 @@ class FF16_Environment : public Environment {
     }
   }
 
-  void set_soil_infiltration_rate(double rate) {
-    soil_infiltration_rate = rate;
-  }
-
   // Core functions
   template <typename Function>
   void compute_environment(Function f_compute_competition, double height_max) {
@@ -92,7 +89,6 @@ class FF16_Environment : public Environment {
   void clear_environment() {
     canopy.clear();
   }
-
 };
 
 inline Rcpp::NumericMatrix get_state(const FF16_Environment environment) {
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 2ce5f7266..5e3ef2683 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -6,11 +6,6 @@
 
 using namespace Rcpp;
 
-#ifdef RCPP_USE_GLOBAL_ROSTREAM
-Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
-Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
-#endif
-
 // Lorenz__ctor
 plant::ode::test::Lorenz Lorenz__ctor(double sigma, double R, double b);
 RcppExport SEXP _plant_Lorenz__ctor(SEXP sigmaSEXP, SEXP RSEXP, SEXP bSEXP) {
@@ -9235,6 +9230,68 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// FF16_Environment__set_extrinsic_driver
+void FF16_Environment__set_extrinsic_driver(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, std::vector<double> x, std::vector<double> y);
+RcppExport SEXP _plant_FF16_Environment__set_extrinsic_driver(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP xSEXP, SEXP ySEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< plant::RcppR6::RcppR6<plant::FF16_Environment> >::type obj_(obj_SEXP);
+    Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP);
+    Rcpp::traits::input_parameter< std::vector<double> >::type x(xSEXP);
+    Rcpp::traits::input_parameter< std::vector<double> >::type y(ySEXP);
+    FF16_Environment__set_extrinsic_driver(obj_, driver_name, x, y);
+    return R_NilValue;
+END_RCPP
+}
+// FF16_Environment__extrinsic_driver_extrapolation
+void FF16_Environment__extrinsic_driver_extrapolation(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, bool extrapolate);
+RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_extrapolation(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP extrapolateSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< plant::RcppR6::RcppR6<plant::FF16_Environment> >::type obj_(obj_SEXP);
+    Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP);
+    Rcpp::traits::input_parameter< bool >::type extrapolate(extrapolateSEXP);
+    FF16_Environment__extrinsic_driver_extrapolation(obj_, driver_name, extrapolate);
+    return R_NilValue;
+END_RCPP
+}
+// FF16_Environment__extrinsic_driver_evaluate
+double FF16_Environment__extrinsic_driver_evaluate(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, double u);
+RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_evaluate(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP uSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< plant::RcppR6::RcppR6<plant::FF16_Environment> >::type obj_(obj_SEXP);
+    Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP);
+    Rcpp::traits::input_parameter< double >::type u(uSEXP);
+    rcpp_result_gen = Rcpp::wrap(FF16_Environment__extrinsic_driver_evaluate(obj_, driver_name, u));
+    return rcpp_result_gen;
+END_RCPP
+}
+// FF16_Environment__extrinsic_driver_evaluate_range
+std::vector<double> FF16_Environment__extrinsic_driver_evaluate_range(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, std::vector<double> u);
+RcppExport SEXP _plant_FF16_Environment__extrinsic_driver_evaluate_range(SEXP obj_SEXP, SEXP driver_nameSEXP, SEXP uSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< plant::RcppR6::RcppR6<plant::FF16_Environment> >::type obj_(obj_SEXP);
+    Rcpp::traits::input_parameter< std::string >::type driver_name(driver_nameSEXP);
+    Rcpp::traits::input_parameter< std::vector<double> >::type u(uSEXP);
+    rcpp_result_gen = Rcpp::wrap(FF16_Environment__extrinsic_driver_evaluate_range(obj_, driver_name, u));
+    return rcpp_result_gen;
+END_RCPP
+}
+// FF16_Environment__get_extrinsic_driver_names
+std::vector<std::string> FF16_Environment__get_extrinsic_driver_names(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_);
+RcppExport SEXP _plant_FF16_Environment__get_extrinsic_driver_names(SEXP obj_SEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< plant::RcppR6::RcppR6<plant::FF16_Environment> >::type obj_(obj_SEXP);
+    rcpp_result_gen = Rcpp::wrap(FF16_Environment__get_extrinsic_driver_names(obj_));
+    return rcpp_result_gen;
+END_RCPP
+}
 // FF16_Environment__canopy_openness
 double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, double height);
 RcppExport SEXP _plant_FF16_Environment__canopy_openness(SEXP obj_SEXP, SEXP heightSEXP) {
@@ -9269,17 +9326,6 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
-// FF16_Environment__set_soil_infiltration_rate
-void FF16_Environment__set_soil_infiltration_rate(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, double rate);
-RcppExport SEXP _plant_FF16_Environment__set_soil_infiltration_rate(SEXP obj_SEXP, SEXP rateSEXP) {
-BEGIN_RCPP
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< plant::RcppR6::RcppR6<plant::FF16_Environment> >::type obj_(obj_SEXP);
-    Rcpp::traits::input_parameter< double >::type rate(rateSEXP);
-    FF16_Environment__set_soil_infiltration_rate(obj_, rate);
-    return R_NilValue;
-END_RCPP
-}
 // FF16_Environment__set_soil_water_state
 void FF16_Environment__set_soil_water_state(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::vector<double> state);
 RcppExport SEXP _plant_FF16_Environment__set_soil_water_state(SEXP obj_SEXP, SEXP stateSEXP) {
@@ -10578,10 +10624,14 @@ static const R_CallMethodDef CallEntries[] = {
     {"_plant_Weibull_Disturbance_Regime__icdf", (DL_FUNC) &_plant_Weibull_Disturbance_Regime__icdf, 2},
     {"_plant_FF16_Strategy__ctor", (DL_FUNC) &_plant_FF16_Strategy__ctor, 0},
     {"_plant_FF16_Environment__ctor", (DL_FUNC) &_plant_FF16_Environment__ctor, 2},
+    {"_plant_FF16_Environment__set_extrinsic_driver", (DL_FUNC) &_plant_FF16_Environment__set_extrinsic_driver, 4},
+    {"_plant_FF16_Environment__extrinsic_driver_extrapolation", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_extrapolation, 3},
+    {"_plant_FF16_Environment__extrinsic_driver_evaluate", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate, 3},
+    {"_plant_FF16_Environment__extrinsic_driver_evaluate_range", (DL_FUNC) &_plant_FF16_Environment__extrinsic_driver_evaluate_range, 3},
+    {"_plant_FF16_Environment__get_extrinsic_driver_names", (DL_FUNC) &_plant_FF16_Environment__get_extrinsic_driver_names, 1},
     {"_plant_FF16_Environment__canopy_openness", (DL_FUNC) &_plant_FF16_Environment__canopy_openness, 2},
     {"_plant_FF16_Environment__clear", (DL_FUNC) &_plant_FF16_Environment__clear, 1},
     {"_plant_FF16_Environment__set_fixed_environment", (DL_FUNC) &_plant_FF16_Environment__set_fixed_environment, 3},
-    {"_plant_FF16_Environment__set_soil_infiltration_rate", (DL_FUNC) &_plant_FF16_Environment__set_soil_infiltration_rate, 2},
     {"_plant_FF16_Environment__set_soil_water_state", (DL_FUNC) &_plant_FF16_Environment__set_soil_water_state, 2},
     {"_plant_FF16_Environment__compute_rates", (DL_FUNC) &_plant_FF16_Environment__compute_rates, 1},
     {"_plant_FF16_Environment__time__get", (DL_FUNC) &_plant_FF16_Environment__time__get, 1},
diff --git a/src/RcppR6.cpp b/src/RcppR6.cpp
index 41566f8ec..672bea1c0 100644
--- a/src/RcppR6.cpp
+++ b/src/RcppR6.cpp
@@ -3752,6 +3752,26 @@ plant::FF16_Environment FF16_Environment__ctor(bool canopy_rescale_usually, int
   return plant::FF16_Environment(canopy_rescale_usually, soil_number_of_depths);
 }
 // [[Rcpp::export]]
+void FF16_Environment__set_extrinsic_driver(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, std::vector<double> x, std::vector<double> y) {
+  obj_->set_extrinsic_driver(driver_name, x, y);
+}
+// [[Rcpp::export]]
+void FF16_Environment__extrinsic_driver_extrapolation(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, bool extrapolate) {
+  obj_->extrinsic_driver_extrapolation(driver_name, extrapolate);
+}
+// [[Rcpp::export]]
+double FF16_Environment__extrinsic_driver_evaluate(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, double u) {
+  return obj_->extrinsic_driver_evaluate(driver_name, u);
+}
+// [[Rcpp::export]]
+std::vector<double> FF16_Environment__extrinsic_driver_evaluate_range(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::string driver_name, std::vector<double> u) {
+  return obj_->extrinsic_driver_evaluate_range(driver_name, u);
+}
+// [[Rcpp::export]]
+std::vector<std::string> FF16_Environment__get_extrinsic_driver_names(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_) {
+  return obj_->get_extrinsic_driver_names();
+}
+// [[Rcpp::export]]
 double FF16_Environment__canopy_openness(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, double height) {
   return obj_->canopy_openness(height);
 }
@@ -3764,10 +3784,6 @@ void FF16_Environment__set_fixed_environment(plant::RcppR6::RcppR6<plant::FF16_E
   obj_->set_fixed_environment(value, height_max);
 }
 // [[Rcpp::export]]
-void FF16_Environment__set_soil_infiltration_rate(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, double rate) {
-  obj_->set_soil_infiltration_rate(rate);
-}
-// [[Rcpp::export]]
 void FF16_Environment__set_soil_water_state(plant::RcppR6::RcppR6<plant::FF16_Environment> obj_, std::vector<double> state) {
   obj_->set_soil_water_state(state);
 }
diff --git a/src/interpolator.cpp b/src/interpolator.cpp
index 6695a27cc..19e76931b 100644
--- a/src/interpolator.cpp
+++ b/src/interpolator.cpp
@@ -18,9 +18,12 @@ void Interpolator::init(const std::vector<double>& x_,
   initialise();
 }
 
-// Compute the interpolated function from the points contained in 'x'
-// and 'y'.
+// Compute the interpolated function from the points contained in 'x' and 'y'.
 void Interpolator::initialise() {
+  // https://stackoverflow.com/questions/17769114/stdis-sorted-and-strictly-less-comparison
+  if (not std::is_sorted(x.begin(), x.end(), std::less_equal<double>())) {
+    util::stop("spline control points must be unique and in ascending order");
+  }
   if (x.size() > 0) {
     tk_spline.set_points(x, y);
     active = true;
@@ -34,6 +37,14 @@ void Interpolator::add_point(double xi, double yi) {
   y.push_back(yi);
 }
 
+// adds point in sorted position (slower than above)
+void Interpolator::add_point_sorted(double xi, double yi) {
+  auto x_upper = std::upper_bound(x.begin(), x.end(), xi); // find smallest number larger than xi
+  x.insert(x_upper, xi); // add xi below that number
+  auto y_upper = std::upper_bound(y.begin(), y.end(), yi);
+  y.insert(y_upper, yi);
+}
+
 // Remove all the contents, being ready to be refilled.
 void Interpolator::clear() {
   x.clear();
@@ -43,6 +54,19 @@ void Interpolator::clear() {
 
 // Compute the value of the interpolated function at point `x=u`
 double Interpolator::eval(double u) const {
+  check_active();
+  if (not extrapolate and (u < min() or u > max())) {
+    auto err = std::string{"Extrapolation disabled and evaluation point of "};
+    err += std::to_string(u);
+    err += " outside of interpolated domain: ";
+    err += "[" + std::to_string(min()) +  ", " + std::to_string(min()) + "]";
+    util::stop(err);
+  }
+  return tk_spline(u);
+}
+
+// faster version of above
+double Interpolator::operator()(double u) const {
   return tk_spline(u);
 }
 
@@ -62,6 +86,10 @@ double Interpolator::max() const {
   return size() > 0 ? x.back() : R_NegInf;
 }
 
+void Interpolator::setExtrapolate(bool e) {
+  extrapolate = e;
+}
+
 std::vector<double> Interpolator::get_x() const {
   return x;
 }
@@ -82,10 +110,10 @@ SEXP Interpolator::r_get_xy() const {
 // points `x=u`, returning a vector of the same length.
 std::vector<double> Interpolator::r_eval(std::vector<double> u) const {
   check_active();
-  const size_t n = u.size();
-  std::vector<double> ret(n);
-  for (size_t i = 0; i < n; ++i) {
-    ret[i] = eval(u[i]);
+  auto ret = std::vector<double>();
+  ret.reserve(u.size()); // fast to do this once rather than multiple times with push_back
+  for (auto const& x : u) {
+    ret.push_back(eval(x));
   }
   return ret;
 }
diff --git a/tests/testthat/test-environment.R b/tests/testthat/test-environment.R
index 4e68dfef4..d1410a716 100644
--- a/tests/testthat/test-environment.R
+++ b/tests/testthat/test-environment.R
@@ -41,3 +41,35 @@ for (x in names(strategy_types)) {
     expect_identical(sapply(hmid, env$canopy$canopy_interpolator$eval), sapply(hmid, interplator$eval))
   })
 }
+
+test_that("FF16 rainfall spline", {
+  context(sprintf("Rainfall-FF16",x))
+  env <- make_environment("FF16w")
+  # get list of extrinsic drivers for the environment
+  expect_equal(env$get_extrinsic_driver_names(), c("rainfall"))
+  
+  # test extrapolation on default spline of y = 1
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 100), 1)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 10000000), 1)
+  
+  # test extrapolation on spline of y = 5.613432
+  env <- make_environment("FF16w", rainfall=5.613432)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 100), 5.613432)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 10000000), 5.613432)
+  
+  ## simple quadratic
+  x <- seq(-10, 10, 0.41)
+  y <- x^2
+  env$set_extrinsic_driver("rainfall", x, y) # overwrites previously created spline
+  
+  # interpolated points
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 2), 4)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", -2), 4)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 3), 9, tolerance=1e-7)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", -3), 9, tolerance=1e-7)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", 5.5), 30.25, tolerance=1e-7)
+  expect_equal(env$extrinsic_driver_evaluate("rainfall", -5.5), 30.25, tolerance=1e-7)
+  
+  ## interpolated range of points
+  expect_equal(env$extrinsic_driver_evaluate_range("rainfall", c(-7, 1, 7.8345)), c(49, 1, 61.37939025), tolerance=1e-6)
+})
diff --git a/tests/testthat/test-interpolator.R b/tests/testthat/test-interpolator.R
index 0a7a587ab..4e621645f 100644
--- a/tests/testthat/test-interpolator.R
+++ b/tests/testthat/test-interpolator.R
@@ -26,6 +26,10 @@ test_that("Splines require sensible data", {
   s <- Interpolator()
   expect_error(s$init(c(1, 2, 3), 1), "Incorrect length input")
   expect_error(s$init(numeric(0), numeric(0)), "insufficient number of points")
+  # ascending order
+  expect_error(s$init(c(3, 2, 4), c(1, 1, 1)), "spline control points must be unique and in ascending order")
+  # unique
+  expect_error(s$init(c(1, 1, 2), c(1, 1, 1)), "spline control points must be unique and in ascending order")
 })
 
 test_that("Spline contains correct data", {
diff --git a/tests/testthat/test-patch.R b/tests/testthat/test-patch.R
index ab3b3379d..e50d7d489 100644
--- a/tests/testthat/test-patch.R
+++ b/tests/testthat/test-patch.R
@@ -42,8 +42,8 @@ for (x in names(strategy_types)) {
     expect_equal(patch$ode_size, env_size)
     
     # either 0 or numeric(0)
-    expect_identical(patch$ode_state, numeric(env_size))
-    expect_identical(patch$ode_rates, numeric(env_size))
+    expect_identical(patch$ode_state, env_state)
+    expect_identical(patch$ode_rates, env_rates)
     
     ## Empty environment:
     patch$compute_environment()
diff --git a/tests/testthat/test-strategy-ff16w.R b/tests/testthat/test-strategy-ff16w.R
index 59c428ec1..bd6cd8d90 100644
--- a/tests/testthat/test-strategy-ff16w.R
+++ b/tests/testthat/test-strategy-ff16w.R
@@ -14,13 +14,15 @@ test_that("FF16w Environment", {
   expect_equal(e$soil$states, 0)
   
   e$compute_rates()
-  expect_equal(e$soil$rates, 0)
+  expect_equal(e$soil$rates, 1) # default rainfall is now y = 1
   
   # Make it rain
-  e$set_soil_infiltration_rate(10)
+  x <- seq(0, 9, 1)
+  y <- rep(5, 10)
+  e$set_extrinsic_driver("rainfall", x, y)
   expect_equal(e$soil$states, 0)
   e$compute_rates()
-  expect_equal(e$soil$rates, 10)
+  expect_equal(e$soil$rates, 5)
   
   # Water logged
   e$set_soil_water_state(100)
@@ -29,7 +31,7 @@ test_that("FF16w Environment", {
   # Check construction
   e <- make_environment("FF16w", 
                         soil_number_of_depths = 1,
-                        soil_infiltration_rate = 10)
+                        rainfall = 10)
   
   expect_equal(e$soil$states, 0)
   e$compute_rates()
@@ -52,31 +54,43 @@ test_that("FF16w Environment", {
   
 })
 
-test_that("Basic run", {
+test_that("Rainfall spline basic run", {
   # one species
   p0 <- scm_base_parameters("FF16w")
   p1 <- expand_parameters(trait_matrix(0.0825, "lma"), p0, FF16w_hyperpar,FALSE)
   p1$birth_rate <- 20
-
+  
   env <- make_environment("FF16w", 
                           soil_number_of_depths = 10,
-                          soil_initial_state = rep(1, 10),
-                          soil_infiltration_rate = 1)
+                          soil_initial_state = rep(1, 10))
+  
+  # init rainfall spline for env
+  x <- seq(0, 110, 0.1)
+  integrand <- function(x) {x^2}
+  y <- integrand(x)
+  env$set_extrinsic_driver("rainfall", x, y)
   
   ctrl <- scm_base_control()
   
   out <- run_scm(p1, env, ctrl)
   
   expect_equal(out$patch$environment$ode_size, 10)
-
+  
+  # model is currently 1/z * int(0, t)(f(x)), where z is depth, 
+  # x is time, f is spline, t is number of years
+  
+  # check the rates are correct, ie 105^2/depth
   expect_equal(out$patch$environment$soil$rates, 
-               c(1.00, 0.50, 0.33, 0.25, 0.20, 0.17, 0.14, 0.13, 0.11, 0.10),
-               tolerance = .01)
+               sapply(seq(1, 10), function(x) { (out$time^2)/x }),
+               tolerance = 1e-7)
   
+  # check the states are correct
+  # model is currently 1/z * int(0, t)(f(x)), where z is depth, 
+  # x is time, f is spline, t is number of years
   expect_equal(out$patch$environment$soil$states,
-               c(105, 52, 35, 26, 21, 17, 15, 13, 11, 10),
-               tolerance = 0.1)
+               sapply(seq(1, 10), function(i) { (1/i) * integrate(integrand, 0, out$time)$value }),
+               tolerance = 1e-5)
   
   expect_equal(out$offspring_production, 16.88946, tolerance=1e-5)
-  expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-5)
-})
+  expect_equal(out$ode_times[c(10, 100)], c(0.000070, 4.216055), tolerance=1e-7)
+})
\ No newline at end of file