model.js

/**** Start of imports. If edited, may not auto-convert in the playground. ****/
var ponds = ee.FeatureCollection("users/kelmarkert/public/ferloPonds"),
    mk_pond = /* color: #d63000 */ee.Feature(
        ee.Geometry.Polygon(
            [[[103.11428283212649, 16.181837537330324],
              [103.11576341150271, 16.181837537330324],
              [103.11509822366702, 16.183238846689054],
              [103.11411117074954, 16.182929735185496]]]),
        {
          "system:index": "0"
        }),
    studyArea =
    /* color: #d63000 */
    /* shown: false */
    ee.Geometry.Polygon(
        [[[-15.866, 16.49],
          [-15.866, 14.193],
          [-12.99, 14.193],
          [-12.99, 16.49]]]),
    chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY"),
    volumne_pt = /* color: #d63000 */ee.Geometry.MultiPoint(),
    lc8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA");
/***** End of imports. If edited, may not auto-convert in the playground. *****/


// Original author: K. Markert
// Based on Soti et al. (2010) --> https://hess.copernicus.org/articles/14/1449/2010/hess-14-1449-2010.pdf
// Edited by: Biplov Bhandari . SCO (4/25/2023)
// This script uses historical CHIRPS for 2019 to derive the Area-Height Relationship.

var elv_org = ee.ImageCollection("projects/servir-wa/SETSM_dem/SETSM_dem2").mosaic().select(["b1"], ["elevation"]);
var elv = elv_org.reproject(ee.Projection('EPSG:4326').atScale(2));
var demScale = elv.projection().nominalScale();
print("DEM resolution", demScale);

Export.image.toDrive({
  image: elv_org,
  description: '2m_DEM',
  folder: '2m_DEM',
  scale: 2,
  crs: 'EPSG:4326',
  maxPixels: 1E13,
});

Map.addLayer(ponds, {}, 'ponds');
Map.addLayer(elv_org, {min: 40, max: 70}, 'dem_2m', false);
// var studyArea = ee.Geometry.Rectangle([-180,-60,180,85])//mk_pond.buffer(10000,100).geometry()

var forecastDays = 365;
// var initDate = ee.Date('2023-12-01');//ee.Date(date);
var initDate = ee.Date('2023-01-01');

var pondId = 75; /// test case was #1

var pond = ee.Feature(ponds.filter(ee.Filter.eq('uniqID', pondId)).first());
print("Sample pond geometry", pond.geometry());
// Map.centerObject(pond, 14);

var rainStudyPeriod = chirps.filterDate(initDate.advance(-20, 'day'),initDate.advance(1, 'day'));
// print('rainStudyPeriod', rainStudyPeriod);
Map.addLayer(ee.Image(rainStudyPeriod.mean()).randomVisualizer(), {}, 'rainStudyPeriod', false);

var initImg = ee.Image(lc8.filterBounds(pond.geometry()).filterDate(initDate.advance(-20, 'day'),initDate.advance(1, 'day')).sort('system:time_start',false).first());
print('initImg', initImg);
var t = ee.Date(initImg.get('system:time_start'));
print('initImg time', t);

// MNDWI on L8 with B3 (green) and B7 (SWIR2)
var initWater = initImg.normalizedDifference(['B3', 'B7']).gt(-0.2); // >-0.2 is a low threshold for MNDWI, potential overestimate of water?
Map.addLayer(initImg.normalizedDifference(['B3', 'B7']), {bands: 'nd', min: -0.4, max: -0.01}, 'MNDWI', false);

var initPct = ee.Number(initWater.reduceRegion({
  geometry: pond.geometry(),
  reducer: ee.Reducer.mean(),
  scale: 30,
  maxPixels: 1e9
}).get('nd'));

Map.addLayer(initImg.select(['B4','B3','B2']),{min:0,max:0.3,gamma:1.3},'Natural-color',true);
Map.addLayer(initImg.normalizedDifference(['B3','B7']).rename("MNDWI"),{min:-0.4,max:-0.1},'Initial Conditions',false);

print("initPct", initPct);

/* ----- model parameterization ----- */
// water balance parameters [from Soti el al. (2010), Table 2]
var k = ee.Image(0.9);            // dimensionless | coefficient expressing soil moisture decrease in time | ranges: 0-1
var Gmax = ee.Image(0.01487);     // m/day | rainfall threshold value start runoff in dry soils | ranges: 0.01-0.02
var L = ee.Image(0.00114);// ***  // m/day | water loss per day | range: 0.005-0.02
var Kr = ee.Image(0.4946);        // dimensionless | runoff coefficient | range = 0.15-0.40
var n = ee.Image(19.89);          // dimensionless | # times catchment area of small pond is larger than the max pond surface area | range: 1-20
var Ac = n.multiply(pond.area()); // sq m | catchement area | range: 0-150,000,000
var alpha = ee.Image(2.514);      // dimensionless | water body shape factor | 1-3


// volume-area-height state variables
var ho = ee.Image(1); // m | pond water height | h0 state --> 1m water height

var _pondMin = elv.reduceRegion({
  geometry: pond.geometry(),
  reducer: ee.Reducer.min(),
  scale:demScale
}).get('elevation');

var pondMin = ee.Image.constant(_pondMin);

print("pond Min height is ", _pondMin.getInfo() + " meters");


var SoInit = ee.Image(0).where(elv.gte(pondMin).and(elv.lte(pondMin.add(ho))), 1);
Map.addLayer(SoInit, {min:0, max:1}, 'SoInit', false);

var nPixels = SoInit.reduceRegion({
  geometry: pond.geometry(),
  reducer: ee.Reducer.sum(),
  scale: demScale
}).get('constant');
print('nPixels', nPixels);

var So = ee.Image(ee.Number(nPixels)).multiply(ee.Image.pixelArea());
Map.addLayer(So, {}, 'So', false);

var SoArea = So.reduceRegion({
  geometry: pond.geometry(),
  reducer: ee.Reducer.first(),
  scale: demScale
}).get('constant');

print('The SoArea is ', SoArea.getInfo() + ' m2');

// calculate initial conditions
var A = ee.Image(pond.area()).multiply(ee.Image(initPct));               // A = pond area at time t [written as A(t)] --> this equation gives the surface area of pond actually covered in water from RS data at given time
var hInit = ho.multiply((A.divide(So)).pow(ee.Image(1).divide(alpha)));       // Eq 6 solve for h(t)
var Vo = (So.multiply(ho)).divide(alpha.add(1));                         // Eq 7 find Vo, Vo is the volume for ho=1m of water height in pond
var vInit = ee.Image(Vo.multiply((hInit.divide(ho)).pow(alpha.add(1)))); // Eq 7 find vInit based on value of Vo just calculated

var currentExtent = A.reduceRegion({geometry:pond.geometry(), reducer:ee.Reducer.first(), scale:30});
print("Maxmum extent of pond:", pond.area()," m2 and current surface area extent: ", currentExtent.get('constant'), " m2");

// set contants
// this was original scale; however this same scale works for CHIRPS coverting value from mm/day to m/day.
var precipScale = ee.Image(1).divide(ee.Image(1e3));

/* ----- start proccessing ----- */

chirps = chirps.select(['precipitation'], ['precip']);
// var precipData = chirps.filterDate(t, t.advance(1, 'day')).filterBounds(studyArea);
// print("Precipitation Data", precipData);

// var dailyPrecip = accumChirps(precipData, t, forecastDays);
var dailyPrecip = accumChirps(chirps, t, forecastDays, precipScale);
print("Daily Precip", dailyPrecip);
Map.addLayer(dailyPrecip, {}, 'dailyPrecip', false);


var pastDays = 30;
// InitIap is a weighted summation of past daily precip amounts used to indicate amount of water in soil
var initIap = calcInitIapWithChirps(chirps.filterDate(t.advance(-pastDays, 'day'), t), pastDays);
// print("Initial Iap", initIap);
Map.addLayer(initIap, {min:0, max:0.1}, "initIap", false);

// set initial conditions with t-1 forcing | this puts all variables into different bands of 1 image for use in volume model
var first = ee.Image(chirps.filterDate(t.advance(-1, 'day'), t).first())
              .multiply(precipScale).addBands(initIap) // converts to m/day units (but is this correct due to s^-1 of original units?)
              .addBands(vInit).addBands(A).addBands(hInit)
              .rename(['precip','Iap','vol','area','height'])//.clip(studyArea)
              .set('system:time_start', t.advance(-1, 'day').millis(), 'system:time_end', t.advance(-1, 'day').millis()).float();
print("initial variables with (t-1) forcing", first);

///////////////////////////////////////////////////////////////////////////////////////////////  // Initialize volume model
var modelOut = ee.ImageCollection.fromImages(dailyPrecip.iterate(accumVolume, ee.List([first])));
print("Model out", modelOut); // why is res of output the res of GFS (~27km2), rather than 30m Landsat res??

Map.addLayer(modelOut.select('vol'), {}, 'modelOut', false);

// returns pond surface area as percentage of total area from pond shapefiles
var pondPct = modelOut.select('area').map(function(img) {
  var pct = img.divide(ee.Image(pond.area())).copyProperties(img, ['system:time_start']);
  return pct;//.where(pct.gt(1),1) --> sets maximum fill at 100%
});
// print("PondPct", pondPct);

print('Temporal Trend of the Volume');
var volTimeSeries = ui.Chart.image.seriesByRegion({
  imageCollection: modelOut.limit(30),
  regions: pond.geometry(),
  reducer: ee.Reducer.mean(),
  band: 'vol',
  scale: demScale,
  xProperty: 'system:time_start',
  seriesProperty: 'label'
});
volTimeSeries.setChartType('ScatterChart');
volTimeSeries.setOptions({
  title: 'Pond volume vs time',
  vAxis: {
    title: 'Volume (m^3)',
  },
  lineWidth: 2,
  pointSize: 3,
  series: {
    0: 'red'
  }
});

print(volTimeSeries);

// export forecastDays
// make image collection before saving
var modelOutLists = modelOut.toList(modelOut.size());
// for (var i=0; i<modelOut.size().getInfo(); i++) {
for (var i=301; i<366; i++) {
  var img = ee.Image(modelOutLists.get(i));
  Export.image.toAsset({
    image: img,
    description: 'img_'+i,
    assetId: 'projects/servir-wa/services/ephemeral_water_ferlo/hydro_model_output/2023/pond_id_' + pondId + '_doy_' + i,
    region: pond.geometry().bounds(),
    scale: demScale,
    maxPixels: 1E13
  });
}

/*---------------------------------------------------------------------------------------*/
// Functions

function accumVolume(img,list) {
  // extract out forcing and state variables
  // "past" equivalent to the x(t-1) state of variables
  var past = ee.Image(ee.List(list).get(-1));//.clip(studyArea);
  var pastIt = past.select('Iap');
  var pastPr = past.select('precip');
  var pastAr = past.select('area');
  var pastHt = past.select('height');
  var pastVl = past.select('vol');
  var nowPr = img.select('precip');//.clip(studyArea);
  var date = ee.Date(img.get('system:time_start'));

  // change in volume model
  var deltaIt = pastIt.add(pastPr).multiply(k);                                // Eq 5
  var Gt = Gmax.subtract(deltaIt);                                             // Eq 4
  Gt = Gt.where(Gt.lt(0),0);                                                   // Eq 4 (cont)
  var Pe = nowPr.subtract(Gt);                                                 // Eq 3
  Pe = Pe.where(Pe.lt(0),0);                                                   // Eq 3 (cont)
  var Qin = Kr.multiply(Pe).multiply(Ac);                                      // Eq 2
  var dV = nowPr.multiply(pond.area()).add(Qin).subtract(L.multiply(pastAr));  // Eq 1 (Qout is assumed to be 0 in Ferlo use case)

  // convert dV to actual volume (add change in volume to the initial volume to get volume at given t step)
  var volume = pastVl.add(dV).rename('vol');
  volume = volume.where(volume.lt(0), 0);

  // empirical model for volume to area/height relationship
  var ht = ho.multiply(volume.divide(Vo).pow(ee.Image(1).divide(alpha.add(1)))).rename('height');
  ht = ht.where(ht.lt(0),0);
  var area = So.multiply(ht.divide(ho).pow(alpha)).rename('area'); //Eq 6
  area = area.where(area.lt(0),1); // constrain area to real values

  // set state variables to output model step
  var step = nowPr.addBands(deltaIt).addBands(volume).addBands(area).addBands(ht)
              .set('system:time_start', date.advance(1, 'day').millis());

  return ee.List(list).add(step.float());
}

// Convert to daily precip
function accumGFS(collection,startDate,nDays) {
  if (nDays>16){
    alert('Max forecast days is 16, only producing forecast for 16 days...');
    nDays = 16;
  }
  var cnt = 1;
  var imgList = [];
  for (var i=0; i<=nDays; i++) {
    var cntMax =(24*(i+1));
    var forecastMeta = [];
    for(cnt;cnt<=cntMax;cnt++){forecastMeta.push(cnt)}
    var dayPrecip = collection.filter(ee.Filter.inList('forecast_hours', forecastMeta));
    imgList.push(dayPrecip.sum().multiply(precipScale)
      .set('system:time_start',startDate.advance(i,'day')));
  }
  return ee.ImageCollection(imgList);
}

function accumChirps (collection, startDate, nDays, scale) {
  if (scale === undefined) { scale = 0 }
  // print("scale accumchirps", scale);
  // chirps has daily values in it so
  var dailyPrecip = ee.ImageCollection(collection.filterDate(startDate, startDate.advance(nDays, 'day')));
  dailyPrecip = dailyPrecip.map(function (img) {
    var sd = img.get('system:time_start');
    var ed = img.get('system:time_end');
    img = ee.Image(ee.Algorithms.If(
      scale !==0,
      img.multiply(scale).copyProperties(img).set('system:time_start', sd, 'system:time_end', ed),
      img.copyProperties(img).set('system:time_start', sd, 'system:time_end', ed)
    ));
    // return img.multiply(precipScale).copyProperties(img).set('system:time_start', sd, 'system:time_end', ed);
    return img;
  });
  return dailyPrecip;
}

function timeScale (img){
  return img.multiply(60*60*6);
}

function accumCFS(collection,s,nDays) {
  var imgList = [];
  for (var i=0; i<nDays; i++) {
    var newDate = s.advance(i,'day');
    var dayPrecip = collection.filterDate(newDate,newDate.advance(24,'hour'))
      .map(timeScale);
    imgList.push(dayPrecip.sum().multiply(precipScale)
      .set('system:time_start',s.advance(i,'day')));
  }
  return ee.ImageCollection(imgList);
}

function calcInitIap(collection, pastDays) {
  var off = pastDays*-1;
  var s = t.advance(off,'day');
  var e = s.advance(pastDays,'day');
  var prevPrecip = collection.filterDate(s,e); // these lines give you precip of select past days (ie past 7 days)

  var dailyPrev = accumCFS(prevPrecip,s,pastDays);

  var imgList = dailyPrev.toList(pastDays);
  var outList = [];  //209-220 Eq 5

  for (var i=0; i<pastDays; i++) {
    var pr = ee.Image(imgList.get(i));
    var antecedent = pr.multiply(ee.Image(1).divide(pastDays-i));
    outList.push(antecedent);
  }
  var Iap = ee.ImageCollection(outList).sum().rename('Iap');
  return Iap;
}

function calcInitIapWithChirps(collection, pastDays) {
  var off = pastDays*-1;
  var s = t.advance(off, 'day');
  var e = s.advance(pastDays, 'day');
  var prevPrecip = collection.filterDate(s, e); // these lines give you precip of select past days

  var dailyPrev = accumChirps(prevPrecip, s, pastDays, precipScale);

  var imgList = dailyPrev.toList(pastDays);
  var outList = [];

  for (var i=0; i<pastDays; i++) {
    var pr = ee.Image(imgList.get(i));
    var antecedent = pr.multiply(ee.Image(1).divide(pastDays-i));
    outList.push(antecedent);
  }
  var Iap = ee.ImageCollection(outList).sum().rename('Iap');
  return Iap;
}


Map.addLayer(pond, {color: 'red'}, 'pond');
// Map.centerObject(pond, 13);