A Semi-Analytical Optical Remote Sensing Model to Estimate Suspended Sediment and Dissolved Organic Carbon in Tropical Coastal Waters Inﬂuenced by Peatland-Draining River Discharges off Sarawak, Borneo

: Coastal water quality degradation is a global challenge. Marine pollution due to suspended sediments and dissolved matter impacts water colour, biogeochemistry, benthic habitats and eventually human populations that depend on marine resources. In Sarawak (Malaysian Borneo), peatland-draining river discharges containing suspended sediments and dissolved organic carbon inﬂuence coastal water quality at multiple locations along the coast. Optical remote sensing is an effective tool to monitor coastal waters over large areas and across remote geographic locations. However, the lack of regional optical measurements and inversion models limits the use of remote sensing observations for water quality monitoring in Sarawak. To overcome this limitation, we have (1) compiled a regional spectral optical library for Sarawak coastal waters, (2) developed a new semi-analytical remote sensing model to estimate suspended sediment and dissolved organic carbon in coastal waters, and (3) demonstrated the application of our remote sensing inversion model on satellite data over Sarawak. Bio-optical data analysis revealed that there is a clear spatial variability in the inherent optical properties of particulate and dissolved matter in Sarawak. Our optical inversion model coupled with the Sarawak spectral optical library performed well in retrieving suspended sediment (bias = 3% and MAE = 5%) and dissolved organic carbon (bias = 3% and MAE = 8%) concentrations. Demonstration products using MODIS Aqua data clearly showed the inﬂuence of large rivers such as the Rajang and Lupar in discharging suspended sediments and dissolved organic carbon into coastal waters. The bio-optical parameterisation, optical model, and remote sensing inversion approach detailed here can now help improve monitoring and management of coastal water quality in Sarawak. SJ extrapolating SIOP datasets bio-optically different regions in between the eastern and western Finally, when all the available SIOPs the retrieved values of TSS exhibited no bias and only a very


Introduction
Knowledge of the distribution of total suspended solids (TSS) and dissolved organic carbon (DOC) is essential for our understanding of water quality degradation and carbon cycle dynamics in coastal oceans. Coastal water quality degradation is a global challenge [1][2][3], which has particularly severe impacts on benthic habitats such as coral reefs [4,5] and seagrass meadows [6,7], and can change the structure, composition and extent of such ecosystems [8][9][10]. This can have knock-on effects on fish populations [11,12],

Field Survey Details
Three fieldtrips were conducted during 2017 to measure bio-optical properties in coastal waters of Sarawak ( Figure 1). The sampling strategy in all the field trips was to collect measurements that represent turbid riverine water, coastal mixing and clear ocean Land cover in Sarawak state is mainly forested area (~64%), of which peat swamp forest occupies 6.1% (Forest Department, Sarawak, 2015). Deforestation in Sarawak for timber production and agriculture development has increased soil erosion in Sarawak rivers [54]. Along the coast in Sarawak, mangrove forests cover 60% of the total coastline. The mangrove forests in Sarawak have suffered from human expansion (agriculture and aquaculture) and pollution through industrialisation and urbanisation [46]. The primary cause of water quality degradation in Sarawak rivers and coastal waters is due to the release of excessive particulate and dissolved substances as a result of changing land cover and land use practices [54,55]. River discharges in Sarawak are known to impact the coastal biogeochemical oceanography in the South China Sea [56].

Field Survey Details
Three fieldtrips were conducted during 2017 to measure bio-optical properties in coastal waters of Sarawak ( Figure 1). The sampling strategy in all the field trips was to collect measurements that represent turbid riverine water, coastal mixing and clear ocean waters. The first field trip was conducted during 8-13 June 2017 (hereafter denoted as SJ). A total of 26 stations were occupied in the suspended sediment-rich waters off the Rajang and Lupar rivers, where coastal peatlands additionally supply terrestrial DOC. Discrete water samples for bio-optical measurements were taken at 20 stations and apparent optical properties were measured during daylight conditions at 10 stations that covered turbid river, coastal and clear ocean locations.
The second field survey was undertaken during 4-10 September 2017 (hereafter denoted as SS). Measurements from this trip represent the ecosystem conditions just before the beginning of the north-east monsoon. A total of 26 field stations in SS covered similar coastal regions of Sarawak as in SJ. Water samples for spectrophotometric analysis (to measure absorption due to particulate and dissolved matter) were collected at 15 stations, backscattering measurements were taken at 24 stations and radiometric measurements made at 19 stations (daytime stations). During both SJ and SS surveys, all the samples were filtered and preserved directly on board the boat.
The third field survey was undertaken during 13-16 September 2017 (hereafter denoted as SM). This field trip covered 23 stations in the western-most area of Sarawak, where suspended sediment concentrations are lower and the peatland-draining Samunsam river provides a discharge of DOC-rich blackwater. Water samples for spectrophotometric analysis were taken at 10 stations, backscattering profiles at 15 stations and radiometric measurements at 16 stations. All stations during SM were sampled from a small outboardpowered boat, and water samples were stored in dark insulated boxes and filtered back on land each afternoon.

Biogeophysical Measurements
Surface water at each station was collected using a clean plastic bucket. Seawater (50-1000 mL) was filtered through pre-rinsed and pre-weighed glass fibre filters (Whatman GF/F, 0.7 µm, 25 mm) to collect Total Suspended Solids. Filters were cleaned with MilliQ water to remove any salt and stored at −20 • C until analysis. Filters were then dried to a constant weight at 75 • C and re-weighed on a microbalance (1 µg readability). The full protocol used to determine TSS is available in [57].
Water for dissolved organic carbon (DOC) analysis was collected in amber glass bottles, and then filtered through a Whatman Anodisc filter (0.2 µm pore size) that was pre-rinsed with 1 M HCl, Milli-Q water, and sample water. The final filtrate (30 mL) was stored in an amber vial, acidified with 0.1 mL of 50% H 2 SO 4 and stored in a refrigerator at 4 • C. Once in the laboratory, DOC measurements were performed on a Shimadzu TOC-L system with a high-salt combustion kit. Instrument performance was monitored using certified deep-sea water from the Hansell Laboratory, University of Miami (42-45 µmol L −1 ). Our analyses yielded a long-term value for the reference water of 47 ± 2.0 µmol L −1 (mean ± 1 SD, n = 51). Procedural blanks prepared in the field almost all contained <10 µmol L −1 , except Remote Sens. 2021, 13, 99 5 of 31 for those prepared in between blackwater river samples, which contained 13-27 µmol L −1 ; a correction for these procedural blanks was not applied [58].

Optical Measurements 2.4.1. Inherent Optical Properties (IOP)
The spectral light absorption coefficient of natural waters (a T ) is the sum of light absorption due to optically active particulate and dissolved substances in the water: a T (λ) = a W (λ) + a P (λ) + a Y (λ) (1) including contributions from water (a W ), particulate matter (a P ) and CDOM (a Y ). In some studies, particulate matter is further separated into contributions from phytoplankton and non-algal particulate matter [59,60]. In this study, the spectral range of λ = 400-700 nm is considered.
To collect CDOM, seawater was filtered as for DOC analysis, and 30-mL samples were preserved with 0.5 mL of a 10% w/v solution of sodium azide (NaN 3 ) solution and stored in amber vials at +4 • C until analysis. While this preservation method follows previous recommendations for remote sensing algorithm development [57], recent research cautions against using NaN 3 as a preservative for CDOM measurements due to blank effects in the ultraviolet range [61]. CDOM samples were analysed within a month of collection. They were brought to room temperature and their absorbance was measured from 230 to 900 nm against a DI water reference using a Thermo Evolution 300 dual-beam spectrophotometer. Instrument performance was checked according to [62]. Reagent blanks of NaN 3 in DI water showed no absorbance over the wavelength range considered here. CDOM spectra were baseline corrected [63] and converted to Napierian absorption coefficients following Equation (2): where a(λ) and A(λ) are, respectively, the absorption coefficient and the absorbance at wavelength λ, and l is the cuvette pathlength in metres. The spectral slope of CDOM (S Y ) was determined by applying a non-linear exponential model to fit a Y between 400 and 700 nm. Spectral CDOM absorption was then modelled as: To measure aP, particles were filtered onto Whatman GF/F filters and stored flat in tissue embedding cassettes in liquid nitrogen until analysis (within three months). The optical density (OD) spectrum was measured over the 350-800 nm spectral range in 0.9 nm increments, using a Cintra 404 UV/VIS dual beam spectrophotometer equipped with an integrating sphere. The OD scans were converted to absorption spectra by first normalising the scans to zero at 750 nm and then correcting for the pathlength amplification using the coefficients of [62]. Unlike CDOM absorption, a P does not have a functional form due to phytoplankton absorption peaks. Particulate absorption has contributions from both non-algal matter and phytoplankton. When dominated by non-algal matter, particulate absorption spectra show an exponential decrease [64], which could be modelled similar to CDOM absorption. However, with increasing phytoplankton presence, the total particulate absorption spectrum will deviate from this exponential behaviour and will include absorption peaks due to phytoplankton pigments (for example, at 443 and 676 nm). To represent such a complex behaviour, in this study spectral particulate absorption a P (λ = 400-700 nm) was therefore modelled as a product of TSS concentration and the mass-specific particulate absorption spectra.
Specific inherent optical properties (SIOP) are computed by normalising IOP to their biogeochemical concentrations. This is useful in understanding and comparing the optical efficiency of particulate and dissolved substances [60]. In this study, mass-specific absorption values were computed for particulate (a P *; m 2 TSS g −1 ) and dissolved substances (a Y * Remote Sens. 2021, 13, 99 6 of 31 (m −1 µM C −1 ). Using these SIOP values for modelling purposes, total absorption could be rewritten as: a T (λ) = a W (λ) + (TSS × a P *(λ)) +(DOC × a Y *(440) × exp(−S Y × (λ−440))) (4) Total light backscattering in the water column can be expressed as: where b bW is the backscattering due to seawater [65] and b bP is the backscattering due to particulate matter suspended in the water column. We measured light backscattering in the water column using a WET Labs ECO BB-9 instrument at wavelengths 412, 440, 488, 510, 532, 595,650, 676 and 715 nm. Raw data from the BB-9 is output as counts from the sensor. By applying the calibration procedure (subtracting the dark offset and multiplying by a scaling factor), the raw counts were converted to meaningful volume scattering coefficients, β T (θ, λ) with units of m −1 sr −1 . Measured β T are further corrected for the absorption in the incident beam following the method suggested for the BB-9 using the spectrophotometric absorption values measured at each station (user guide, WET Labs ECO BB-9, 2017). The corrected volume scattering coefficients comprise contributions from particles and water in the backward direction (θ = 117 • ). Particulate volume scattering (β P ) was obtained by subtracting the volume scattering of water [65]. Particulate backscattering coefficients, b bP (λ) with units of m −1 , were then determined as: For a diversity of water types, Boss and Pegau [66] suggested a value of 1.1 for the factor X, to convert scattering in an angle in the back direction to the backscattering coefficient (full details in WET Labs ECO BB-9, 2017).
The spectral slope of b bP (λ) = TSS × b bP *(650) × (λ/650) SbbP (7) In this study, total backscattering was then modelled as: Backscattering albedo (uT) is the ratio of total backscattering to the sum of total backscattering and total absorption [67,68]. Backscattering albedo directly influences the remote sensing reflectance and the light availability in the water column [69][70][71]. We computed the backscattering albedo using the BB-9 measured backscattering data and our spectrophotometric measurements of total absorption as:

Apparent Optical Property Measurements
In this study, hyperspectral (320-950 nm at 2 nm resolution) radiometric measurements were made using Trios-RAMSES radiometers (https://www.trios.de/en/ramses.html). Apparent optical properties (AOP) of the water column such as above-surface downwelling irradiance Es(0+, λ), underwater downwelling irradiance Ed(depth, λ) and upwelling radiance Lu(depth, λ) were measured. Deployment of instruments followed the protocols suggested in [72]; the radiance sensor had a field of view of 7 • . Measurements were made only during high sun elevation conditions (10 AM to 4 PM during the day) under stable sun and sky conditions.

Relationship between Underwater Remote Sensing Reflectance and Backscattering Albedo
Remote sensing reflectance just below the surface, r rs (0−) is known to vary as a function of single scattering albedo, as described by the following equation [67,71]: where f indicates the relationship between irradiance reflectance and IOPs ( [74] and references therein) and Q (ratio of upwelling radiance to irradiance) [75] represents the directional nature of upwelling light field. The term f/Q is equivalent to G0(uT) + G1(uT), and r rs (0−) could be rewritten as below: r rs (0−) = (G0 × (uT) + G1 × (uT)) × uT (13) where G0 and G1 represent the combined influence of illumination, sea surface conditions and marine volume scattering function [73]. For a given r rs (0−), uT could also be derived by rewriting the above equation as: Measurements of r rs (0−) and uT from this study were used in the equations above to estimate regional G0 and G1 values for Sarawak coastal waters.

Statistics
We describe the variation in bio-optical properties using their minimum and maximum values, and their mean and percent coefficient of variance (% CV = 100 (standard deviation/mean)). Relationships between parameters were assessed using linear correlation (r). Strength of regression was described using the coefficient of determination (R 2 ). The performance of optical models was evaluated using bias and mean absolute error (MAE), which were computed following Seegers et al. [76]. We used log-transformed data when computing bias and MAE, as the data span multiple orders of magnitude. For interpretation purposes, they were converted out of log10 space. Bias values closest to unity are considered least biased and values less than unity indicate a negative bias. MAE indicates the average magnitude of model errors, where the increase above unity indicates the relative measurement error (i.e., MAE = 1.00 means no error, while 1.10 means 10% error). Mean percent relative deviation (MPRD(%) = mean(100 (model − measurement)/measurement) was used to compare the simulated and measured optical signature (uT or R rs ).

Remote Sensing Data
Satellite remote sensing data from the MODIS Aqua sensor were used to demonstrate the application of the remote sensing inversion model. MODIS Aqua has bands at 412, 443, 488, 531, 555, 667, and 678 nm in the visible region. We downloaded Level-2 reflectance products from https://oceancolor.gsfc.nasa.gov/, which are atmospherically corrected [77] and are reported to have good performance in multiple water types [78]. However, in some pixels in coastal waters close to land, the atmospheric correction fails (resulting in negative remote sensing reflectance), which are flagged and removed before the inversion model is applied. For demonstration purposes, relatively cloud-free MODIS Aqua images were used.

Model Structure
The remote sensing inversion model developed in this study has three main components ( Figure 2): (1) a regional spectral optical library of specific inherent optical properties, (2) a forward model to simulate total absorption, total backscattering and backscattering albedo, and (3) an inversion mechanism to derive biogeochemical parameters through spectral matching of remote sensing and forward-modelled backscattering albedo. The model is sufficiently flexible to be applied to both in-situ and satellite-measured remote-sensing reflectance. The model components are described in detail in Sections 3.2-3.6.

Regional Spectral Optical Library
We compiled a regional Sarawak spectral optical library (SSOL; Figure 3) for use in the optical model. The data used for this model was collected during June and September 2017, which represent coastal ecosystem conditions during different monsoon periods. SSOL mainly includes specific inherent optical properties that help simulate total absorption, total backscattering and backscattering albedo as described in Equations (4), (8) and (9).
Our measured aP*(λ) varied over a wide range ( Figure 3A), with higher values in the blue region and relatively lower values in the green-red region. Similar to coastal waters

Regional Spectral Optical Library
We compiled a regional Sarawak spectral optical library (SSOL; Figure 3) for use in the optical model. The data used for this model was collected during June and September 2017, which represent coastal ecosystem conditions during different monsoon periods. SSOL mainly includes specific inherent optical properties that help simulate total absorption, total backscattering and backscattering albedo as described in Equations (4), (8) and (9). vary most strongly spatially rather than temporally in coastal waters of Sarawak. This is consistent with previous analysis of the mixing behaviour of DOC and CDOM, which found that there was a single, strong relationship between DOC and CDOM in our dataset across all seasons [58]. Notably, Martin et al. [58] also included DOC and CDOM data for both the eastern and western area of our study region during March 2017, representing the end of the northeast monsoon (when other optical properties were not measured). The CDOM and DOC dataset used here was also used in [58] to analyse the biogeochemistry of dissolved organic matter in Sarawak, but that study reported other CDOM parameters than we use here (e.g., aY(350) and spectral slopes at UV wavelengths). Light backscattering due to particulate matter also varied spatially more than temporally in Sarawak coastal waters ( Table 1). Values of bbP varied in response to changing TSS concentrations (linear correlation, r = 0.87) in coastal waters, with lower values in clear marine waters and higher values in turbid sediment-rich coastal waters near river mouths. Spatially, the mean bbP(650) values were higher in the eastern region off the Rajang delta, and were slightly higher during September than June, although the range of values largely overlapped between the two seasons (Table 1). Regionally, bbP values measured in Sarawak were comparable to those measured in the Mahakam delta [82], but were higher by an order of magnitude compared to values from the Berau estuary, Indonesian Borneo (bbP(660) = 0.0345) [85], despite the relative proximity of the Mahakam and Berau rivers. The spectral shape of bbP followed a hyperbolic function with spectral slopes (SbbP) varying between 0.40 and 3.75 nm −1 , with mean values that are comparable to other coastal regions [59,85,86]. Higher values of SbbP were found in the more TSS-rich eastern region than in the western region (Table 1). Higher values of SbbP indicate smaller average particle diameters [86], and such small particles should offer larger surface area for photon interaction and increase the particulate backscattering efficiency. However, in this data set, we only observed a poor correlation (r = 0.25) between bbP*(650) and SbbP, suggesting that other particle characteristics (such as shape, composition and porosity; [87]) influenced the particulate backscattering efficiency. The backscattering efficiency covered a large and mostly Our measured a P *(λ) varied over a wide range ( Figure 3A), with higher values in the blue region and relatively lower values in the green-red region. Similar to coastal waters elsewhere [64,68,[79][80][81], the spectral shape of a P * varied due to the changing composition of particulate matter. In clear waters, with mainly phytoplankton as particulate matter, a P * spectra clearly showed absorption peaks around 443 and 676 nm. In contrast, increasing concentrations of non-algal particulates masked the phytoplankton absorption peaks in more turbid coastal waters, yielding a P * spectra that decrease towards the red region. Due to the presence of phytoplankton absorption peaks, the a P *(λ) could not be written in a functional form; for this reason, the full-spectrum a P * values were included in SSOL. The spectral library does not explicitly separate particulate absorption into phytoplankton and non-algal particulate matter, as our focus was on distinguishing between CDOM and particulate absorption. However, the spectral library could easily be expanded to include subgroups of particulate matter if necessary.
Sarawak coastal waters are overall rich in CDOM ( Table 1). The highest CDOM absorption was found in the Samunsam and Sematan river mouths in the western part of study area (measured during SM2017; Figure 1). In the coastal waters influenced by the Rajang River, mean values of a Y (440) were similar in June and in September 2017. Our CDOM absorption values for Sarawak are similar to values reported from the Mahakam Delta (Indonesian Borneo) and surrounding coastal waters (mean a Y (440) = 0.87 m −1 ) [82], but on average higher than the values measured off the east coast of Peninsular Malaysia both during the NE monsoon (0.15 m −1 ) and at other times (0.05 m −1 ) [80]. The spectral slope of CDOM is a good indicator of the origin of CDOM, with higher values (about 0.021 nm −1 ) associated with marine sources and lower values (about 0.014 nm −1 associated with terrestrial sources) [83]. In addition to these two groups, in Sarawak coastal waters we also observed spectral values lower than 0.01 nm −1 in 6% of the samples (occurring both in clear and turbid waters). Analysis of fluorescence spectra has shown that DOC and CDOM in the coastal waters of our study region was to a large degree terrestrial in origin, with peatlands representing the main terrestrial source [84]. The absorption efficiency (as indicated by a Y *) of CDOM varied by two orders of magnitude in the western part of the study area, where both the blackwater Samunsam River and relatively clear marine waters were sampled (Table 1; SM2017). In the eastern part of the study area, a Y * covered a similar range during both the June and September field trips, which was smaller than in the western region. Overall, this dataset indicates that the CDOM absorption properties vary most strongly spatially rather than temporally in coastal waters of Sarawak. This is consistent with previous analysis of the mixing behaviour of DOC and CDOM, which found that there was a single, strong relationship between DOC and CDOM in our dataset across all seasons [58]. Notably, Martin et al. [58] also included DOC and CDOM data for both the eastern and western area of our study region during March 2017, representing the end of the northeast monsoon (when other optical properties were not measured). The CDOM and DOC dataset used here was also used in [58] to analyse the biogeochemistry of dissolved organic matter in Sarawak, but that study reported other CDOM parameters than we use here (e.g., a Y (350) and spectral slopes at UV wavelengths).
Light backscattering due to particulate matter also varied spatially more than temporally in Sarawak coastal waters ( Table 1). Values of b bP varied in response to changing TSS concentrations (linear correlation, r = 0.87) in coastal waters, with lower values in clear marine waters and higher values in turbid sediment-rich coastal waters near river mouths. Spatially, the mean b bP (650) values were higher in the eastern region off the Rajang delta, and were slightly higher during September than June, although the range of values largely overlapped between the two seasons (Table 1). Regionally, b bP values measured in Sarawak were comparable to those measured in the Mahakam delta [82], but were higher by an order of magnitude compared to values from the Berau estuary, Indonesian Borneo (b bP (660) = 0.0345) [85], despite the relative proximity of the Mahakam and Berau rivers. The spectral shape of b bP followed a hyperbolic function with spectral slopes (S bbP ) varying between 0.40 and 3.75 nm −1 , with mean values that are comparable to other coastal regions [59,85,86]. Higher values of S bbP were found in the more TSS-rich eastern region than in the western region (Table 1). Higher values of S bbP indicate smaller average particle diameters [86], and such small particles should offer larger surface area for photon interaction and increase the particulate backscattering efficiency. However, in this data set, we only observed a poor correlation (r = 0.25) between b bP *(650) and S bbP , suggesting that other particle characteristics (such as shape, composition and porosity; [87]) influenced the particulate backscattering efficiency. The backscattering efficiency covered a large and mostly overlapping range between the two regions and seasons, but was on average slightly higher in the western region (SM; Figure 1), and in June compared to September in the eastern region (Table 1). These results indicate that backscattering properties of particulate matter show especially spatial variability in Sarawak coastal waters, and bio-optical models for this region must therefore capture this variability. The high spatial variability in the inherent optical properties of Sarawak reflects partly the diversity of river discharges: rivers originating in peatlands carry high CDOM and DOC concentrations but lower TSS concentrations, while rivers draining mineral soils carry high TSS but lower DOC and CDOM concentrations. It also reflects the strong on-shore to off-shore gradients in TSS and CDOM, as particles settle out beyond river plumes, and CDOM is diluted with water from the open sea and photobleached. These results highlight the limitation of extrapolating optical properties from one region to another. To incorporate this process-and region-specific variability, we categorised the spectral optical library into three optical water types (OWT): clear marine waters, turbid coastal waters, and highly turbid coastal waters (influenced by river discharges). Our classification of OWT was based on the reflectance signature: similar to Aurin et al. [71], we used the reflectance ratio to distinguish between highly turbid "red" waters (Type-3: R rs (650):R rs (443) > 1.5), coastal waters (Type-2: R rs (650):R rs (443) ≤ 1.5), and clear marine waters (Type-1: R rs (650):R rs (443) ≤ 0.05). Once categorised, each set of SIOPs (a P *, a Y * and b bP *) is supplied to the optical model as a set and not as independent values, because they are associated together with an individual process (for example, a river plume). This treatment of sets of SIOPs helps to represent the optical properties in a realistic fashion, avoiding the random mixing of SIOPs that were never encountered in combination. Thus, each set of specific inherent optical properties was categorised and associated with one OWT. Such a categorisation helps the selection of realistic SIOP sets for inverse modelling.

Conversion of Underwater Remote Sensing Reflectance and Backscattering Albedo
Remote sensing reflectance just below the surface, r rs (0−), varies as a function of single scattering albedo [67,71,88,89], and the coefficients G0 and G1 differ between clear ocean waters and turbid coastal waters (see Section 2.4.1). Fixed values of G0 = 0.0949 and G1 = 0.0794 were suggested for clear ocean waters by Gordon et al. [88] and G0 = 0.0895 and G1 = 0.1247 for coastal waters by Lee et al. [67]. In contrast, Aurin et al. [71] proposed a set of spectrally varying f/Q values for conversion (see Table 3 in [71]).
In Sarawak coastal waters, the parameterisations by Gordon et al. [88], Lee et al. [67], and Aurin et al. [71] generally underestimated uT values, even though they showed acceptable linearity, with the spectrally varying conversion factors of Aurin et al. [71], yielding slightly lower bias than the other two approaches ( Figure 4A-C). To improve agreement between modelled and measured uT values, we used Sarawak region-specific values (G0 = 0.0508) and G1 = −0.0271), which were generated through optimisation of regional values of r rs (0−) and uT. Using these Sarawak-specific values in the relationship between r rs (0−) and uT (Equation (14)) resulted in a better match of uT ( Figure 4D).

Forward Modelling of Backscattering Albedo
The forward model simulates total absorption, total backscattering and backscattering albedo spectra based on the SIOP values in the spectral library. CDOM absorption is simulated as the product of DOC-specific CDOM absorption at 440 nm and DOC concentration, and is then combined with the CDOM spectral slope to yield a simulated CDOM spectrum (Equation (3)). Particulate absorption and backscattering are simulated as the products of TSS concentration and spectral TSS-specific particulate absorption, and TSS-specific backscattering (Equations (4) and (8)). Total backscattering spectra are then computed with the particulate backscattering spectral slope and by adding the backscattering of water. Total backscattering and absorption were then used in the forward model to compute the backscattering albedo (uT(λ); Equation (9); Figure 2). This model performed well across all water types that we sampled, with uT spectra showing peaks in the blue region (400-500 nm) for clear marine waters and shifting towards the red region (600-700 nm) in turbid coastal waters ( Figure 5A-D). The range in spectral shapes of modelled uT is overall consistent with observations from eastern Borneo [79,82] and other similar tropical coastal waters, such as Singapore [90] and off the Amazon [91].

Forward Modelling of Backscattering Albedo
The forward model simulates total absorption, total backscattering and backscattering albedo spectra based on the SIOP values in the spectral library. CDOM absorption is simulated as the product of DOC-specific CDOM absorption at 440 nm and DOC concentration, and is then combined with the CDOM spectral slope to yield a simulated CDOM spectrum (Equation (3)). Particulate absorption and backscattering are simulated as the products of TSS concentration and spectral TSS-specific particulate absorption, and TSSspecific backscattering (Equations (4) and (8)). Total backscattering spectra are then computed with the particulate backscattering spectral slope and by adding the backscattering of water. Total backscattering and absorption were then used in the forward model to compute the backscattering albedo (uT(λ); Equation (9); Figure 2).

Inversion Model
Our inversion modelling approach uses the Levenberg-Marquardt (LM) algorithm [93] to iteratively change spectral SIOPs (from SSOL) for the selected OWT, and then to adjust the concentrations of TSS and DOC in the forward model until the best fit is found between the supplied and predicted backscattering albedo (uT(λ)). Briefly, the LM algorithm is an iterative technique that locates the minimum of a multivariate function that is expressed as the sum of squares of non-linear real-valued functions [94]. As we have coupled the forward model (backscattering albedo) with the iterative mechanism of the LM, our inversion model will be referred to as iSAM (iterative semi-analytical model). We used In the eastern region off the Rajang delta, both in June and September 2017, stations in relatively clear marine waters had higher uT values in the blue region (400-500 nm) due to lower concentrations of TSS and DOC ( Figure 5A,B). Closer to the coast, the concentrations of TSS and DOC were higher, which reduced uT in the blue due to increased absorption by particulate and dissolved substances, and increased uT in the red region due to higher backscattering. Stations close to and in the Rajang River distributary mouths had the highest concentrations of DOC and TSS in the eastern region, and uT spectra consequently had the highest values in the red region owing to the higher backscattering values ( Figure 5A,B). The forward-modelled uT in the western region (SM) showed some differences to the eastern region ( Figure 5A,B versus Figure 5C). For example, in clear marine waters (solid blue lines in Figure 5), uT values were on average 22% higher in the western region (SM) than in the eastern region for similar concentrations of TSS and DOC. These higher values of uT in the western region (clear marine waters) are due to higher backscattering efficiency of particles (b bP *(650) = 0.0157) relative to the eastern region in September (b bP *(650) = 0.0028). In contrast, closer to river mouths (solid red line; Figure 5B,C), higher DOC concentrations in the western region reduced the uT values in the blue-green (400-600 nm) part of the spectrum while the presence of TSS resulted in high values in the red (600-700 nm). Moreover, the shape and magnitude of uT spectra close to the river mouths changed according to hydrodynamic conditions ( Figure 5D). Sarawak experiences a meso-tidal range of up to 5 m, which together with the multiple river discharges results in optically complex conditions, as also seen in coastal waters elsewhere [60,92]. In June 2017, we sampled coastal waters off the Igan distributary in the Rajang delta (2.875 • N 111.638 • E) over a period of 14 h to understand the influence of Igan river discharges. Initially at a salinity of 13.1, the values of TSS and DOC were 25.4 mg L −1 and 209.3 µM C, which caused higher uT values in the red region (dashed red line; Figure 5D). As salinity decreased to 8 and TSS concentrations increased by 51%, uT increased across the spectrum (solid red line; Figure 5D). The river plume receded after another 10 h: salinity increased to 21.5 and TSS and DOC were reduced, such that the uT spectrum shifted to have a maximum in the green region (500-600 nm; solid green line in Figure 5D). Because our spectral library includes SIOPs from different hydrodynamic conditions encountered in Sarawak coastal waters, it thus has the ability to simulate the range of optical signatures encountered in these dynamic coastal waters.

Inversion Model
Our inversion modelling approach uses the Levenberg-Marquardt (LM) algorithm [93] to iteratively change spectral SIOPs (from SSOL) for the selected OWT, and then to adjust the concentrations of TSS and DOC in the forward model until the best fit is found between the supplied and predicted backscattering albedo (uT(λ)). Briefly, the LM algorithm is an iterative technique that locates the minimum of a multivariate function that is expressed as the sum of squares of non-linear real-valued functions [94]. As we have coupled the forward model (backscattering albedo) with the iterative mechanism of the LM, our inversion model will be referred to as iSAM (iterative semi-analytical model). We used an IDL implementation of the LM algorithm ('MPFIT' code available at http://purl.com/net/ mpfit; [95]), which we applied to our backscattering albedo model with input data of SIOPs. LM algorithms have been successfully used with in-situ data in clear ocean waters [96], complex coastal waters [97] and also with satellite remote sensing data [98].
To test the iSAM model, we grouped the dataset according to the three field expeditions (SJ, SS and SM), which represent spatially and temporally varying conditions. We first calibrated the model with the spectral library from each data subset individually and applied it to predict TSS and DOC for the remaining two datasets by supplying the computed backscattering albedo from those datasets. We then calibrated the model with the full dataset and applied it to predict the TSS and DOC for the full dataset from the backscattering albedo, to quantify the inherent model error.
For TSS, we found that iSAM did not perform very well when the spectral library from only one data subset was used to predict TSS for the other two subsets, with either relatively large biases, mean absolute errors, or both ( Figure 6A-C). Model performance was relatively better when using the spectral library from either SJ or SS to predict TSS in the other two data subsets ( Figure 6A,B) and was worst when using only the spectral library from SM to predict TSS in SJ and SS ( Figure 6C). These results demonstrate the limitations of extrapolating SIOP datasets between bio-optically different regions in Sarawak coastal waters, especially between the eastern and western region. Finally, when all the available SIOPs were supplied to iSAM, the retrieved values of TSS exhibited no bias and only a very small MAE ( Figure 6D). The small MAE (5%) reflects the inherent error associated with the inversion modelling (iSAM) approach when retrieving TSS values in Sarawak region.  We found similar results for DOC: when applying the spectral library of just one individual data subset (SJ, SS, or SM) to predict DOC for the other two subsets, iSAM performed relatively poorly, showing large biases and MAEs ( Figure 7A-C). The fact that SIOPs from SS (or SJ) could not very well predict DOC concentrations for SJ (or SS) uT spectra might suggest that there were seasonal differences in the extent of biogeochemical transformations of DOC and CDOM in the eastern region, especially from photobleaching of terrestrial CDOM [58]. However, it may well simply reflect the fact that only a limited range of bio-optical conditions were sampled during each trip, with the SJ dataset especially not capturing clearer marine waters well (Table 1). However, when the spectral library contained all available SIOPs, it retrieved concentrations of DOC with good accuracy, showing no bias and only minimal MAE. This indicates only a small inherent model error for DOC retrieval (8%). Similarly, only small errors were noted for intermediate modelling components, such as a Y (440), a P (440) and b bP (650) (see Figure A1A-C). Overall, our results for DOC and TSS highlight the importance of using region-and processes-specific SIOPs in Sarawak to achieve better retrievals of TSS and DOC concentrations.

Predictive Errors of Inversion Model
We estimated the predictive error of the iSAM model using the leave-one-out crossvalidation (LOO-XV) method. LOO-XV is a K-fold cross-validation technique that uses

Predictive Errors of Inversion Model
We estimated the predictive error of the iSAM model using the leave-one-out crossvalidation (LOO-XV) method. LOO-XV is a K-fold cross-validation technique that uses most of the data for calibration [99]. In this study, we used n−1 samples for calibration and used the left-out sample for validation. This method is repeated until all the samples have been used for validation. Measured and predicted values for each station were then used to compute the bias and MAE. These predictive errors represent the error that results if the measured biogeophysical and optical parameters for every station are predicted without using the SIOPs for that station, but instead using the SIOPs from the station with the closest-matching uT value. Since all SIOPs in the spectral library are obviously available to iSAM during normal use, the predictive errors derived from LOO-XV are almost certainly a considerable overestimate of the real errors associated with this model. However, the error estimate from LOO-XV can be taken as indicative of the errors associated with using iSAM in cases where the true SIOPs are not well represented by the existing spectral library. For cases where the true SIOPs are well matched by an SIOP combination in our spectral library, the true error is given by the inherent model error (see Figures 6D and 7D) rather than by the predictive error.
The predictive error of the iSAM model for TSS estimation was low in Sarawak coastal waters ( Figure 8A), when tested with in-situ optical measurements. Inversion model performance varied between subsets, with bias = 0.72 and MAE = 1.63 for SJ, and with bias = 0.76 and MAE = 1.91 for SS. However, in the SM subset the model performance was relatively poor (bias = 2.52 and MAE = 2.57), mainly due to stations in SM region that were influenced by the Samunsam river discharges. At these stations, large values of DOC (382.6, 652.8, 1798.7 µM C) were noted and the resulting average model bias was 10.78 for TSS retrieval. Stations with such high concentrations of DOC were very close to river mouths and adjacent land, and such pixels are generally flagged as problematic in remote sensing imagery anyway. In case such river discharges are transported further out to sea without dilution, then the inversion model may be limited in its performance in such pixels. iSAM performance varied in changing water conditions: in clear marine waters (TSS < 1 mg/L) the model bias was 1.83, in coastal waters (1 < TSS < 100 mg/L) the model bias was 1.06, and in highly turbid waters (TSS > 100 mg/L), iSAM under-estimated TSS values with an average bias of 0.29. Highly turbid waters with TSS concentrations above 100 mg/L were only encountered within rivers and their estuaries, but not in coastal waters beyond river mouths, so this underestimate is unlikely to affect coastal remote sensing applications of iSAM. Overall, the iSAM inversion model has good predictive error statistics in coastal waters, and relatively poor performance in clear marine waters and highly turbid river discharges. The performance of iSAM in Sarawak coastal waters could be improved by including more bio-optical measurements from clear and highly turbid waters in SSOL.
Retrieval of DOC with the iSAM model shows good agreement with DOC measurements ( Figure 8B) with low predictive errors (bias = 0.96, MAE = 1.36). Consistent with the performance on the full dataset, iSAM also performed well in subsets SJ (bias = 1.09, MAE = 1.21), SS (bias = 1.21, MAE = 1.31) and SM (bias = 0.62, MAE = 1.61). The largest deviations were observed in the SM dataset for the high-DOC stations (>300 µM DOC) taken directly adjacent to the river mouth. The mismatches in DOC estimates indicate the bio-optical variability of DOC (and CDOM) present in the coastal waters along the eastern coast of Sarawak [58,84]. Comparison of iSAM performance across different DOC gradients also showed minimal predictive errors in waters with DOC < 100 µM DOC (relatively clear waters; bias = 1.10, MAE = 1.28) and for waters where 100 > DOC < 300 µM C (coastal waters; bias = 1.07 and MAE = 1.20). The highest predictive error (bias = 0.37 and MAE = 2.66) was observed for stations in highly turbid coastal waters at the river mouths. In general, iSAM thus performs well in Sarawak coastal waters for the retrieval of DOC concentrations. However, like TSS retrieval, its performance deviates in highly turbid conditions created by direct river discharges close to shore. Predictive errors associated with the IOPs are shown in Figure A2.   Figure 9B) and red wavelengths ( Figure 9C). Figure 9D demonstrates the influence of short-term variability of SI-OPs on simulated uT at the same geographic location (2.87°N/111.63°E) within just 14 h. The solid line in Figure 9D is from station SJ15 (salinity = 8.0), the dashed line indicates uT simulated using SIOPs from SJ14 (04:50 h prior to SJ15; salinity = 13.0), and the dotted line indicates SIOPs from SJ16 (09:20 h after SJ15; salinity = 21.5). This observed variability in uT indicates the changes in the bio-optical properties at this location in response to hydrodynamic conditions. The examples shown in Figure 9 demonstrate the importance of using region-and process-specific SIOPs to simulate realistic optical signatures (uT). This means that the more SIOPs from different hydrodynamic and biogeochemical conditions in Sarawak are included in SSOL, the better are the chances that iSAM will invert the signature successfully. Currently, SSOL has 38 full spectral SIOP sets. To improve the performance of bio-optical models, future studies should measure and improve the Sarawak spectral optical library.   Figure 9 demonstrate the importance of using region-and process-specific SIOPs to simulate realistic optical signatures (uT). This means that the more SIOPs from different hydrodynamic and biogeochemical conditions in Sarawak are included in SSOL, the better are the chances that iSAM will invert the signature successfully. Currently, SSOL has 38 full spectral SIOP sets. To improve the performance of bio-optical models, future studies should measure and improve the Sarawak spectral optical library. Remote Sens. 2021, 13, x FOR PEER REVIEW 21 of 32

Model Application to Satellite Remote Sensing Data
To demonstrate the use of the iSAM inversion model for remote sensing purposes, we applied it to a MODIS AQUA image from 14 June 2017. Cloud cover is a major challenge for optical remote sensing in equatorial regions, and still affected some areas of this rare, largely cloud-free, image. The true-colour image ( Figure 10A) shows clear blue marine waters in the South China Sea, turbid sediment-rich waters all along the coast, and turbid river discharges near multiple river mouths. Also visible is the dark DOC discharge emanating from the Maludam peatland region.

Model Application to Satellite Remote Sensing Data
To demonstrate the use of the iSAM inversion model for remote sensing purposes, we applied it to a MODIS AQUA image from 14 June 2017. Cloud cover is a major challenge for optical remote sensing in equatorial regions, and still affected some areas of this rare, largely cloud-free, image. The true-colour image ( Figure 10A) shows clear blue marine waters in the South China Sea, turbid sediment-rich waters all along the coast, and turbid river discharges near multiple river mouths. Also visible is the dark DOC discharge emanating from the Maludam peatland region. The MODIS AQUA Level 2 product that we supplied to iSAM was first converted from remote sensing reflectance at the bottom of the atmosphere to underwater remote sensing reflectance (Equation (11)) and further to backscattering albedo (uT(λ); Equation (14)). Cloud-free pixels were classified into the three optical water types (OWT, see Section 3.2) according to the reflectance signature, allowing suitable SIOPs to be selected from the spectral library and supplied to iSAM for inversion. This classification discriminated well between different water types, clearly revealing an eddy feature with river-influenced water (OWT2) moving offshore that is not evident in the true-colour image ( Figure 10A,B). Also observed is the presence of highly turbid waters (OWT3) in narrow regions very close to the coast. Unfortunately, most areas with OWT3 pixels were masked in the Level 2 product due to the failure of the atmospheric correction in these waters. Figure 10C shows the spectral matching achieved by iSAM for this MODIS Aqua image. Mean percent relative deviation (MPRD) is the mean of spectral deviations between MODIS measurements and the iSAM-simulated uT(λ). In general, iSAM achieved good agreement in most of the coastal waters (OWT2), with MPRD varying between ±10%.
Relatively higher levels of model performance were achieved in OWT2 because the SSOL includes sufficient SIOP members that represent the hydrodynamic and biogeochemical conditions in Sarawak coastal waters. Closer to the land, MPRD reached up to −20%, while in clear marine water (OWT1) the model performance varied between −20 to −40%. The larger deviations observed in model performance in OWT1 indicate a lack of sufficient SIOP members representing these waters in our spectral library.
TSS and DOC concentrations derived using iSAM for the image of 14 June 2017 in coastal water types are within the range we measured ( Figure 10C,D; Table 1). TSS values were relatively higher in coastal waters than compared to offshore regions, with highest values seen in pixels closer to river mouths (Rajang, Igan, Lupar, Sematan and Lundu). It is also evident from the image that the eastern coast of Sarawak had higher concentrations of TSS delivered by several rivers (e.g., Rajang and Lupar) than compared to the western coast of Sarawak. The MODIS Aqua-derived DOC concentrations showed similar patterns to TSS, with generally high values throughout the coastal waters and peak values closest to river mouths ( Figure 10E). In offshore marine waters (OWT1), the model clearly under-estimates DOC concentrations, because the spectral library lacks data for CDOM-poor marine DOC that would predominate in waters unaffected by terrestrial DOC inputs. An analysis of fluorescence spectra measured in parallel with our DOC and CDOM data suggests that even in the clearest waters we sampled, 20-40% of DOC was of terrestrial origin [84], and therefore characterised by higher specific absorption than typical for marine waters. An interesting feature to note is the elevated DOC concentration in the eddy structure, which appears to carry DOC from the peatland-rich Rajang delta further into the ocean. Such remote sensing images of the distribution of DOC and TSS can supply important information for our understanding of carbon and sediment transport and processing in complex coastal ocean regions of Sarawak. Figure 11 shows TSS and DOC estimates from MODIS Aqua for different time periods in 2017, revealing the high spatio-temporal variation in DOC and TSS in Sarawak. In these equatorial coastal waters, it was difficult to find completely cloud free images. For this demonstration purpose, we selected relatively cloud-free images from 23 January, 4 May and 21 May (all 2017). Figure 11A,B, from 23 January, show conditions during the north-east monsoon period, with generally high concentrations of TSS and DOC in coastal waters. In particular, strong river discharges carrying TSS and DOC could be observed in front of the Sadong river mouth and all across coastal waters in front of the Rajang and Igan river mouths. However, the concentrations are relatively low in coastal waters along the eastern side (Samunsam-Lundu region) of Sarawak. Figure 11C,D show MODIS images from 4 May during the early south-west monsoon period. In these two images, TSS and DOC concentrations are high all along the Sarawak coast. The higher concentrations in coastal waters could be linked to strong river discharges, especially from the Rajang delta, and the Lupar and Samunsam rivers. The Rajang and Igan along the eastern coast had higher concentrations and a wider spread of discharges compared to the rivers in the western region. Higher DOC concentrations delivered by the Samunsam river are consistent with our in-situ measurements from this blackwater river (Table 1). Similar to the image from January, the western region has the largest areal extent of waters with relatively low concentrations of TSS and DOC, indicating higher water quality. Finally, images from 21 May show higher concentrations in front of the Rajang and Igan river mouths, but a reduced spread of these high concentration features out to sea. Also, the strong river discharge with high concentrations of TSS and DOC from the Samunsam river was not present, probably owing to differences in tidal conditions during the satellite overpass.
Overall, the TSS and DOC images derived using iSAM presented in Figure 11 clearly demonstrate the strength of remote sensing maps in understanding Sarawak coastal waters. These iSAM-derived images reveal hotspots of river discharge that contribute particulate and dissolved substances to the coastal ecosystems. Figure 10E further shows that hydrodynamic processes can transport particulate and dissolved substances far into the ocean. Overall, it appears that the western region has higher water quality in terms of the concentrations of dissolved and particulate terrestrial substances, consistent with the presence of coral reefs in that region. riod, with generally high concentrations of TSS and DOC in coastal waters. In particular, strong river discharges carrying TSS and DOC could be observed in front of the Sadong river mouth and all across coastal waters in front of the Rajang and Igan river mouths. However, the concentrations are relatively low in coastal waters along the eastern side (Samunsam-Lundu region) of Sarawak. Figure 11C,D show MODIS images from 4 May during the early south-west monsoon period. In these two images, TSS and DOC concentrations are high all along the Sarawak coast. The higher concentrations in coastal waters could be linked to strong river discharges, especially from the Rajang delta, and the Lupar and Samunsam rivers. The Rajang and Igan along the eastern coast had higher concentrations and a wider spread of discharges compared to the rivers in the western region. Higher DOC concentrations delivered by the Samunsam river are consistent with our in-situ measurements from this blackwater river (Table 1). Similar to the image from January, the western region has the largest areal extent of waters with relatively low concentrations of TSS and DOC, indicating higher water quality. Finally, images from 21 May show higher concentrations in front of the Rajang and Igan river mouths, but a reduced spread of these high concentration features out to sea. Also, the strong river discharge with high concentrations of TSS and DOC from the Samunsam river was not present, probably owing to differences in tidal conditions during the satellite overpass.
Overall, the TSS and DOC images derived using iSAM presented in Figure 11 clearly demonstrate the strength of remote sensing maps in understanding Sarawak coastal waters. These iSAM-derived images reveal hotspots of river discharge that contribute particulate and dissolved substances to the coastal ecosystems. Figure 10E further shows that hydrodynamic processes can transport particulate and dissolved substances far into the ocean. Overall, it appears that the western region has higher water quality in terms of the concentrations of dissolved and particulate terrestrial substances, consistent with the presence of coral reefs in that region.
2. An optical inversion model coupled with the Sarawak spectral optical library successfully retrieved TSS and DOC concentrations from in-situ measurements of backscattering albedo.

3.
Remote sensing demonstration products of MODIS Aqua-derived TSS and DOC show the influence of river discharges on Sarawak coastal waters.
The bio-optical analysis presented in this study extends our understanding of coastal bio-optical properties and their variability in response to changing hydrodynamic processes. The optical model and its parameterisation will benefit both remote sensing and coastal ecosystem modelling efforts in complex coastal waters such as Sarawak. The remote sensing inversion products developed in this study now provide a synoptic view of processes affecting coastal water quality of Sarawak. New spatial and temporal remote sensing products, as demonstrated in this study, will help improve coastal monitoring practices and management strategies in Sarawak coastal waters. This is the first study in Borneo that uses a comprehensive set of biogeochemical and optical properties. However, there are still further improvements that could help expand our understanding of coastal biogeochemical and optical dynamics in Sarawak. Expanding the Sarawak spectral library further across Borneo (e.g., into Brunei and Sabah) will help large-scale spatial studies along the western coast of Borneo. The bio-optical datasets included in SSOL currently only represent the conditions during June and September 2017, and it might improve the model solution if the SSOL is further expanded with more samples that represent the peak periods of north-east and south-west monsoons. Finally, in this study we validated the inversion model mainly with in-situ optical data. It is a challenge to achieve satellite match-ups in tropical coastal regions due to the extensive cloud cover, but it is important to attempt further validation match-ups in these waters. We therefore recommend that the local Sarawak research community should conduct further bio-optical studies to help improve satellite match-up datasets in this region.