A method for quantifying freshwater discharge rates from satellite observations and Lagrangian numerical modeling of river plumes

Satellite remote sensing is an efficient tool for identifying buoyant river plumes in the coastal ocean that are formed by the interaction between river discharge and ambient seawater. A new method for reconstructing the volume rates of river discharge based on the shape, extent and orientation of plumes is described that combines the output from a Lagrangian numerical model and analyses of satellite imagery. At the first step in the procedure, a high resolution satellite image is used to identify the river plume. The spatial characteristics of the plume as seen in the image are not determined solely by the current river discharge rate, as they also depend on the hydrographic features in the sea and atmospheric forcing. A previously developed and validated hydrodynamic model for river plumes is run with a variety of forcing conditions to identify the discharge rate that provides the best match between modeled and observed plumes. The method can be applied to estimate indirectly discharge from small rivers and streams, many of which lack direct measurements. Here it has been applied and validated against in situ data for two rivers feeding the eastern part of the Black Sea.


Introduction
Continental discharge is a significant component of the global hydrological cycle, although it accounts for no more than one tenth of the input part of the ocean water budget on the global average (Oki and Kanae 2006). In many instances, the buoyant fluvial inflow causes surface density stratification over large shelf areas, and plays an important role in physical, chemical, and biological processes (Garvine 1987, Mallin et al 2005, Boyer et al 2006, Korotenko 2007. In particular, the high nutrient content of river discharge together with the stratification induced by plumes in winter and early spring may lead to exceptionally early algal blooms and may cause eutrophication and hypoxia over broad shelf areas (Wiseman and Kelly 1994, Rabailais 2010, Tang et al 2003, Zhou et al 2008. Even small river discharges often form significant buoyant plumes in the areas adjacent to river estuaries (Devlin et al 2001, Zhurbas et al 2011. River plumes, and especially those generated by small rivers, are characterized by high temporal variability caused by external forcing (Yankovsky and Chapman 1997, Fong and Geyer 2002, Gan et al 2008 and mixing processes (Simpson andBowers 1984, Nash 1994). River plumes are usually bounded by distinct frontal zones (Zavialov et al 2003, Horner-Devine et al 2015. Terrigenic constituents such as river-borne suspended matter, which are readily observed using airborne and satellite sensors, can act as natural 'markers' of plumes (Walker 1996, Klemas 2012, 2013. The spatial scales and patterns, as well as the dynamic behavior of a plume, hold information about the incoming river discharge and can be detected in satellite imagery to provide quantitative insights into discharge volume. Data on river discharge rates are essential for many practical purposes. Nowadays, however, less than 60% of continental discharges at the global scale are subject to regular measurements at the river mouths (Fekete et al 1999). Moreover, the number of gauge stations and the availability of discharge data have been decreasing constantly since the 1980s , The Ad Hoc Group 2001. Regular direct measurements of discharge are performed only for a limited number of rivers, generally major ones or those flowing through populated, easily accessible areas. Most of the numerous small rivers with catchment areas below 10 000 km 2 do not have any gauge stations, even though their basins cover about 20% of the area draining into the oceans, and their total contribution to the influx of the global continental discharge is estimated at about 25% Syvitski 1992, Fekete et al 2002). Small rivers are also believed to bring about 40% of the total sediment load to the ocean (Milliman et al 1999). At the regional scale, their contribution can be much more significant.
In some cases, the missing in situ discharge data can be substituted by remote airborne and satellite measurements (Bjerklie et al 2003, Syed et al 2010. Several indirect methods of estimating river discharge based on hydrological modeling and remote measurements by radar altimeters (Papa et al 2010, Tarpanelli et al 2013, interferometric synthetic aperture radars (Smith andPavelsky 2008, Hirpa et al 2013), and passive microwave radiometers (Brakenridge et al 2012) have been developed since the 1980s. Indirect methods based on remote sensing hold promise of partly substituting laborious and expensive direct discharge measurements (Barrett 1998). However, conversion of data collected by remote measurements into river discharge values requires additional atmospheric or hydrological information, such as river bed morphology, bathymetry (Birkett et al 2002, Alsdorf et al 2007 or land-atmosphere water mass fluxes (Syed et al 2009), which is not always easily available. Besides, most of these methods have been applied only to relatively large rivers, whilst small-size rivers with narrow river beds remain largely unaddressed, due to both the limited spatial resolution of satellite sensors (Ducet et al 2000, Njoku et al 2003, Enjolras et al 2006 and the lack of necessary hydrological information about small rivers.
In this article, we present a method for estimating river discharge that combines satellite imagery and numerical modeling. The novelty of this method comes from a new approach based on remote sensing and subsequent modeling of a buoyant plume formed by the river discharge in the sea, in contrast to estimating the volume of the river flow itself as proposed in most previous studies. River plumes tend to diffuse over a wide area adjacent to the river estuary, unlike the river stream bounded by the banks. This fact is crucial for satellite observations of small-size rivers. Accordingly, the method presented here can be applied to small rivers with narrow riverbeds, as well as complex watershed systems like deltas and coastal marshes. In particular, our algorithm was adapted and validated for the small-size rivers of the eastern part of the Black Sea region.
The article is organized as follows. Section 2 provides the general description of the model. Section 3 is focused on adaptation and validation of the model including detailed information about the study region and the remote sensing data used. The discussion and the conclusions are presented in section 4.

Method
The sequence of the method is shown diagrammatically in figure 1. Firstly, the spatial extent and the shape of the plume generated by the river discharge under consideration is identified using satellite imagery of the coastal zone adjacent to the river estuary. Secondly, the inverse problem of obtaining the discharge rate of the river (or rivers) that forms the examined plume is solved. For this purpose, a series of numerical simulations of the river runoff spread are performed under a variety of prescribed external forcing conditions (e.g., wind, coastal circulation, tides, waves) and river discharge parameters (e.g., river inflow velocity, river discharge rate). We assume that the results of the river plume simulation performed under the real forcing conditions (including the sought-for real discharge), i.e., the 'modeled' plume, will be similar to the plume observed in satellite image, i.e., the 'satellite' plume. By varying the forcing parameters in the model, we evaluate the accordance between the model results and the satellite image to achieve the best match. As a result, we obtain the estimate of the river discharge as one of the parameters of the best-fitting configuration. The model consists of three sequential procedures, their description is given below.

Detection of river plume borders
The first procedure identifies the borders of a buoyant plume in the input satellite image ( figure 1(a)). Buoyant river plumes and surrounding coastal waters possess very different physical and chemical properties. Therefore different ocean surface characteristics measured by satellite instruments can be employed to distinguish buoyant plumes in coastal areas.
Previous related studies focused mainly on elevated or reduced concentrations of different constituents of water in river plumes that can be identified by optical satellite measurements. These substances can be regarded as passive tracers of freshwater discharge if the temporal scale characteristic for their sedimentation, transformation or consumption by biota is much longer than the temporal scale of the plume dissipation (James 2002). Surface concentrations of chlorophyll a (Chl-a), total suspended solids (TSS), and colored dissolved organic matter (CDOM) can form sharp and persistent gradients at the plume borders visible in satellite imagery at high or medium spatial resolution (e.g., 30 m for OLI/Landsat-8, 100 m for OLCI/Sentinel-3, 250 m for MERIS/Envi-Sat). A number of previous studies successfully mapped river plumes including plumes formed by relatively small rivers (Nezlin and Digiacomo 2005, Schroeder et  The most common way to define the plume border is prescribing a threshold value for a certain water characteristic. The corresponding isoline can then be regarded as the borderline (Walker 1996). However, concentrations of tracers in both river and sea waters are subject to significant synoptic and seasonal variability. Therefore, the threshold values that define isolines distinguishing plume and sea water have to account for this variability and be chosen individually for different images of the same river plume depending on season and other details of image acquisition. This difficulty can be partly circumvented by considering values of gradient of the quantity used as the indicator of the plume border (e.g., the line of maximum gradient). Another approach implemented in the present model consists in cluster analysis of the coastal area using several surface water characteristics as plume indicators (Oliver et al 2004, Thomas and Weatherbee 2006, Palacios et al 2009, Bainbridge et al 2012. Among the clusters identified through this algorithm, those adjacent to the river mouth will indicate the area occupied by the river plume.

Numerical simulation of river plume spreading
The second procedure of the model simulates the spreading of river discharge under prescribed settings received as input data ( figure 1(b)). These data consist of the static parameters of the model domain (shoreline geometry, location and geometry of the estuary, Coriolis parameter, river and sea water salinity, etc) and dynamical external forcing conditions (wind speed and direction, background coastal current velocity, river discharge). The ultimate objective of the model is to find dynamical parameters that provide the best agreement between the 'modeled' plume and the 'real' plume seen in the satellite image. Even relatively small variations of wind stress, coastal circulation, and river discharge rate can significantly change plume behavior (Kingsford and Suthers 1994, Devlin et al 2001, Gaston et al 2006, Zhurbas et al 2011. Therefore, to ensure an accurate search for parameters that provide agreement between the 'modeled' and 'satellite' plumes, we need to perform a large number of model runs, varying the parameters of external forcing with small increments for every single case. This would imply a high computational cost of the model and limit applicability of the method. This problem was solved in the following way. Firstly, instead of widely used Eulerian models, we employ here a Lagrangian model called STRiPE (Surface-Trapped River Plume Evolution), recently developed specifically for simulating river plumes (Osadchiev and Zavialov 2013). This model tracks the motion and dissipation of imaginary individual 'particles' of the river water discharged into the sea. The main advantage of STRiPE lies in its ability to provide realistic results at relatively low computational cost as compared to Eulerian models. The general description of the STRiPE model and some results of its implementation for simulating river plumes in the Black Sea and the South China Sea can be found in Osadchiev and Zavialov (2013) and Korotenko et al (2014).
Secondly, for every processed image of the river plume, we perform a large number of model runs under different model configurations, sequentially changing all external forcing parameters with relatively small increments. The ranges and increments of the external forcing parameters used in the simulations are presented in table 1. The resulting set of 'modeled' plumes is organized as a 'catalog', storing both the contour of a plume and its model configuration. The creation of the catalog with a large number of 'modeled' plumes under various forcing conditions (6×8×5×2×16=7680 model experiments, see table 1) was possible due to the computational efficiency of STRiPE. This approach provides a marked advantage when analyzing multiple satellite images of the plume taken on different days. Once the database is constructed, the procedure of identifying the bestfitting discharge rate for every individual satellite scene reduces to a simple search for the best match in the catalog without having to perform any additional model simulations.
In most cases, the database contained an element with sufficient level of agreement with the considered 'satellite' plume, which was then used as the first-order approximation for a final iterative procedure. Because the identified forcing parameters already provide relatively good agreement between the 'satellite' and 'modeled' plumes, for 'fine tuning' of the discharge estimate we can vary the values of the input parameters within only narrow ranges around this firstorder approximation. In addition, we can assume that the river plume area is a monotonic function of the input parameters within these narrow ranges (which is not true otherwise). Hence, we can use a bisectional iterative procedure for the input parameters in the prescribed ranges to estimate river discharge more accurately.

Comparison of 'satellite' and 'modeled' plumes
The third procedure of the model compares the 'satellite' plume border calculated by the first procedure and the 'modeled' plume border simulated by the second procedure ( figure 1(c)). This procedure is used to identify the best-fitting element in the database and in the final iterative procedure. For this purpose, we need to define a numerical function that measures the level of similarity between two representations of the plume. On the one hand, we model salinity distribution produced by interaction of sea and river waters, in which case we consider the salinity-derived 'modeled' plume. On the other hand, we identify plume area in satellite imagery and deal with concentrations of TSS, CDOM and Chl-a and thereby obtain the 'satellite' plume. Fortunately, in most cases (although not always) plumes defined through salinity and the tracer constituents practically coincide. The third procedure is based on two functions, which compare plume areas and plume shapes. The first function simply calculates the difference in the plume areas. The second function evaluates similarity of the contours of the two considered plumes. This function calculates scalar products of contour-vectors of plumes, thus evaluating similarity between the contours of two considered plumes. A detailed description of the second function is given in the appendix.

Validation
The presented method was applied and tuned for the Mzymta River (annual average discharge 49 m 3 s −1 ) and the Sochi River (15 m 3 s −1 ) situated in the eastern part of the Black Sea coast. Some background information about these rivers can be found in Jaoshvili (2002) (figure 2).
Gauge measurements at these rivers generally are performed with daily frequency, while the other rivers of the study region are not covered by any regular measurements of discharge. Validation runs against in situ data were performed for 77 different days during a period of one year (April 2011-April 2012) for the Mzymta River and for 53 days during a period of 7 months (April 2011-October 2012) for the Sochi River. The validation period for the Sochi River is 1 The angle is counted from the local direction of the shoreline at the river mouth. Owing to the small sizes of the considered river plumes and proximity of the shore, we assume only alongshore currents; therefore there are only two variants of current direction. 2 One relative unit is equal to estimated mean monthly discharge of the simulated river.
shorter than for the Mzymta River because the gauge data of the Sochi River during winter and early spring in 2010/2011 are not available. As the input remote sensing data, we used satellitederived maps of TSS, CDOM and Chl-a retrieved from MERIS/EnviSat satellite imagery (300 m spatial resolution) provided by the European Space Agency. The data were processed using the BEAM Toolbox software (Fomferra and Brockmann 2014). Surface concentrations of TSS, CDOM, and Chl-a were retrieved from Level 1 satellite products using MERIS Case-2 Regional water processing module based on neural network algorithms Schiller 2007, 2008).
The first procedure of the model that identifies borders of river plumes performs cluster analysis of the coastal area using satellite-retrieved concentrations of TSS, CDOM and Chl-a as the input data. The outcomes from this procedure are contours of the areas adjacent to the river estuaries along the sea shore united in the 'plume' cluster ( figure 3).
The choice of clustering algorithm, metrics for distance, and the number of clusters significantly influence the clustering results and therefore has to be based on a strong theoretical background. This is especially true for the choice of a clustering algorithm that depends on properties of the analyzed data sets. In this work we used an objective method (k-means clustering), that is less biased as it is based on patterns naturally occurring in the data, as opposed to clustering algorithms based on neural networks or fuzzy logic (Chang et al 2002, Oliver et al 2004. The choice of a centroid-based clustering method was due to the fact that the river plume area in the analyzed satellite data is represented by a separate convex set as far as it has significantly different origin than the surrounding ocean waters. The satellite data sets were standardized to a [0; 1] range, and the distances between the elements were calculated using an Euclidian metric to reproduce the homogeneity of the dimensions of the clustering space. The number of clusters and the number of iterations were chosen to provide stable results.
The k-means clustering algorithm described well the typical patterns of the river plumes and showed satisfactory agreement with the in situ data collected in the areas adjacent to the mouths of 12 different rivers in the eastern part of the Black Sea during ten monitoring surveys of the Shirshov Institute of Oceanology in -2012.
The model domain covered the coastal regions adjacent to the Mzymta and Sochi river estuaries correspondingly that are influenced by the river plumes. Spatial resolution of the model was prescribed as 25 m, while the model integration time step was set to 10 min. The density of ambient sea waters was set to 1017 kg m −3 , while the river water density was equal to 1000 kg m −3 . The Coriolis parameter for the study region was selected as 1.028×10 −4 s −1 corresponding to the latitude 43.5°N. The external forcing conditions were varied as in the table 1. The Mzymta and Sochi river plumes are characterized by relatively small sizes (less than 10 km) and are formed at the narrow and steep shelf. Coastal circulation at the upper layer of the study area was therefore regarded as alongshore currents. The maximal wind speed forcing for the numerical simulations was set equal to 5 ms −1 because 93% of 10-minutes averaged wind speeds registered at the meterological station at the study region did not exceed 5 ms −1 . Simulation ranges and increments of the varied river discharge depended on the mean seasonal river size to optimize the number of model runs. For this purpose river discharge rates for numerical experiments were prescribed in relative units instead of m 3 s −1 .
Seventy seven validation days during a 12-month period (April 2011-April 2012) for the Mzymta River and 53 validation days during a 7-month period (April 2010-October 2010) for the Sochi River were selected according to availability of cloud-free satellite data (figure 4). Discharge for both rivers is characterized by significant seasonal and synoptic variability, the latter is caused by active precipitation events at the study region. Both validation selections cover all seasons, although summer season data prevail for the Mzymta River, while the data for the Sochi River are more uniformly distributed within the year.
The indirect estimates of the discharge rate for the Mzymta River and Sochi River showed good agreement with the in situ measurements performed at gauge stations. They reproduce the synoptic and seasonal variability of river discharge reasonably well (figures 4 and 5). The overall comparison between the satellite-derived discharge and gauge measurements for the 130 trial cases is illustrated by figure 5. The mean absolute error of 15% quantifies the ability of the model to reproduce absolute values of measured discharge. The root mean square error is 26.9 m 3 s −1 caused by the sensitivity of this characteristic to large errors. We note that the model showed better agreement with gauge measurements for low discharge cases compared to flooding conditions. This feature is presumably due to the following reason: the value of river discharge is fixed within the individual simulation of a river plume spread, that is, the method outputs an integrated value of discharge rate for a short period of time preceding the time shown at the analyzed satellite image. As a result, it shows worse agreement with daily gauge measurements for high discharge periods that are mainly induced by active, but short-term precipitation events that result in high variability of the discharge rate ( figure 4). The coefficient of determination R 2 is 0.95, which shows that the   model has good level of accuracy in reproducing the intra-annual variability of river discharges.

Discussion and conclusions
A new method of estimating river discharge using satellite imagery and numerical modeling was developed. The general idea of this method is a conversion of satellite-derived properties of river plumes in the sea into information about the river discharge. The advantage of this approach compared to other satellite methods reported earlier lies in the fact that most commonly used satellite sensors cannot resolve smallsize rivers, but can resolve the buoyant plumes formed by these rivers. The proposed method can be used to evaluate freshwater discharges from small rivers with narrow riverbeds and complex watershed systems like deltas and coastal marshes where other indirect methods encounter substantial difficulties.
The presented method strongly depends on the availability and quality of satellite imagery. Firstly, it cannot be used during ice coverage periods as river plumes are formed beneath sea ice, hindering their detection by remote sensing. Secondly, coastal areas obscured by clouds often accompany rainfall-induced flooding events especially for rivers with small watershed basins and a quick response of discharge to precipitation events. This fact limits the applicability of optical satellite imagery for river runoff estimation. However, this problem could be overcome using nonoptical satellite products which do not depend on cloud coverage, for instance SAR products, which also have a very good spatial resolution. The significantly different surface dynamics of river plumes and ambient ocean results in different surface roughness within plume and sea waters, which can then be used to detect river plumes as was shown in Zheng et al (2004), DiGiacomo et al (2004, and Jiang et al (2009). Finally, the availability of satellite imagery at the area under consideration determines the frequency of runoff estimations (except during ice coverage periods). Once satellite images of the proper quality are received on a daily basis, this method is able to reproduce annual variations of river discharge. Otherwise, sparse, but regular estimations of river runoff provided by the method and averaged on the seasonal or annual period can be used to evaluate inter-seasonal or inter-annual discharge variability.
Small and medium-size rivers are responsible for a considerable part of the overall continental runoff into the ocean. They significantly influence the everyday life of millions of people living at the sea coast through their impact on water quality, pollution, fishery, and engineering activities (Emmett et al 2006, Milliman et al 2007, Syvitski and Saito 2007. Thus, remote sensing techniques aimed at obtaining indirect estimates of the discharge from small and medium-size rivers hold promise for providing improved quantitative assessments and new insights into the land-ocean fluxes of freshwater, dissolved and suspended matter, nutrients and pollutants.
Future development of this method will be aimed at extending the described technique to estimate the total discharge of ungauged rivers and streams at the sea or regional scale. Nowadays, more than 35% of the world's population experiences water shortage, and about 80% are threatened by an insecure water supply (Gleick 2003, Vörösmarty et al 2010. Therefore, the presented method when extended to regional scales could be extremely valuable for many part of the world with limited or absent direct gauge measurements caused by the remoteness of the river locations and/or poor economic and social national conditions. In order to compare 'satellite' and 'modeled' plumes we need to use a function that evaluates the similarity of two plume contours. This function deals with the contour-vectors of plumes, i.e., vectors of complex numbers, derived as follows. Once we have a continuous sequence of pixels that form the border of a plume in a satellite image, we can set a starting point at this curve (e.g., the river mouth) and a direction of motion along the contour (e.g., clockwise). Then we start to move along the curve and for every step towards a new contour pixel we can assign a complex number a+bi, where a and b stand for the length of the step in zonal (x) and meridional (y) directions, respectively. As far as we consider unit steps, there are only eight variants for a+bi, namely, 1, −1, i, −i, 1+i, −1+i, 1 −i, −1−i. We can define the scalar product q of two contour-vectors α=(α 1 , K, α n ) and β=(β 1 , K, β n ), if they have the same length n (the number of elements), in the following way: q , , = + + -The scalar product of two vectors is a product of the lengths of these vectors and a cosine of an angle between them. The large (small) angle between the vectors corresponds to the small (large) value of their scalar product; it can therefore be used as a measure of similarity between two vectors. For the same reason, the scalar product of two contour-vectors can be used as a measure of similarity between the two contours. We describe the normalized scalar product r of two contour-vectors in the following way: r , where , i n i i 1 ( ) a a a = å = is a norm of a contourvector.
Due to the Cauchy-Schwarz inequality for any normalized scalar product r, 0r1 and r=1 if and only if α=kβ, where k is a complex number. The module of the product of two complex numbers is equal to the product of their modules, and the argument of the product is equal to the sum of the arguments. This means that contour kβ is a rotated and scaled contour β, the magnitudes of rotation and scaling are determined by the number k. As a result, the normalized scalar product of contour-vectors resized to one length can be employed as a numerical criterion of similarity between contours of two plumes. More detailed information about the applied mathematical tools is given in Han and Fan (1994).