Analysis of Water Yield Changes in the Johor River Basin, Peninsular Malaysia Using Remote Sensing Satellite Imagery

Abstract

Changes in land-use–land-cover (LULC) affect the water balance of a region by influencing the water yield (WY) along with variations in rainfall and evapotranspiration (ET). Remote sensing satellite imagery offers a comprehensive spatiotemporal distribution of LULC to analyse changes in WY over a large area. Hence, this study mapped and analyse successive changes in LULC and WY between 2000 and 2015 in the Johor River Basin (JRB) by specifically comparing satellite-based and in-situ-derived WY and characterising changes in WY in relation to LULC change magnitudes within watersheds. The WY was calculated using the water balance equation, which determines the WY from the equilibrium of precipitation minus ET. The precipitation and ET information were derived from the Tropical Rainfall Measuring Mission (TRMM) and moderate-resolution imaging spectroradiometer (MODIS) satellite data, respectively. The LULC maps were extracted from Landsat-Enhanced Thematic Mapper Plus (ETM+) and Landsat Operational Land Imager (OLI). The results demonstrate a good agreement between satellite-based derived quantities and in situ measurements, with an average bias of ±20.04 mm and ±43 mm for precipitation and ET, respectively. LULC changes between 2000 and 2015 indicated an increase in agriculture land other than oil palm to 11.07%, reduction in forest to 32.15%, increase in oil palm to 11.88%, and increase in urban land to 9.82%, resulting in an increase of 15.76% WY. The finding can serve as a critical initiative for satellite-based WY and LULC changes to achieve targets 6.1 and 6.2 of the United Nations Sustainable Development Goal (UNSDG) 6.

1. Introduction

Satellite-based spatio-temporal land-use land-cover (LULC) change analysis can provide comprehensive and widespread water yield (WY) distribution information [1]. However, changes in land use affect the water balance of a region by influencing WY along with variations in rainfall and evapotranspiration (ET) [2]. Land use patterns can directly change ecosystem types, landscape patterns, and ecological processes [3]. These have a significant influence on the WY on the surface, a crucial parameter for planning and managing raw water supply [4,5]. This affects ecosystem services such as biodiversity, WY characteristics, and rainfall [6].

The conventional WY measurement is based on a plot representation within a catchment/watershed. The plot could be observations of a single or a few points to measure rainfall and ET for a given period, producing a WY value based on the representation of the point or the entire polygon. However, with advances in remote sensing (RS) technologies, WY information could be obtained spatially and temporally at the pixel level through specific processing. The advantages of using RS to extract WY information include (1) comprehensiveness for every pixel; (2) reduced costs—minimum field observations are made for calibration and validation only; (3) feasibility even for remote and ungagged watersheds; and (4) derivation of WY from satellite using selected hydrologically established methods, i.e., the water balance equation, which provides a simple and robust approach. However, the approach is replicable at any location, provided there is available corresponding local satellite data and the needed processing parameters [7].

Previous studies established that satellite-based precipitation data offer high accuracy worldwide [8,9]. The higher-level data products of the Tropical Rainfall Measuring Mission (TRMM) and Global Precipitation Measurement (GPM) data products of the satellite-based precipitation data provide daily rainfall estimation within a three-hour interval at a high spatial resolution of 0.25°, which can be localized to reach a higher resolution [10,11]. Several other studies have raised concerns about the effects of changes in land cover due to LULC changes in the WY [12,13]. An assessment of the relationship between land use and hydrological cycle components was also carried out [14,15,16,17]. Runoff modelling for drought and flood monitoring was reported [18]. These studies, as mentioned above, only widely reported the effects of LULC changes on WY, without taking cognizance of the influence of individual LULC classes on WY. This could lead to an inadequate understanding of WY dynamics.

Between 1984 and 2015, about 90,000 km² of permanent surface water was lost from the Earth’s surface, primarily due to human extraction and increased evapotranspiration caused by climate change [19]. Thus, accounting for water is critical to developing policies and procedures that ensure the maintenance of a regional or basin-scale water balance over time [19,20]. Although Malaysia is considered a low-water-stress country based on the World Resource Institute Ranking [21], Johor, despite its good river system [22] centred around the Johor River, which provides 60% of its water supply [23], is known for water scarcity [24]. According to the National Water Resource study, demand for water in Johor has consistently increased, and this trend will continue [25]. This may not be unconnected to the increasing population, urbanisation, industrialisation and large-scale agricultural activities [26], increasing demands for water [27].

The increasing human activities in Johor are also causing large-scale LUCCs [23,28,29], which, in turn, are affecting the water balance with intermittent and regular incidences of droughts and flooding [27,29,30,31,32,33]. At present, the city-state of Singapore draws 40% of its water from the Johor River [34]. Hence, the Johor River Basin (JRB) plays a pivotal role in the water resource management of both Peninsular Malaysia and the island state of Singapore. Consequently, this study achieved precise satellite-based reporting of changes in WY from six LULC classes (agriculture, bare land, forest, oil palm, urban, and water bodies) with specific objectives: (a) to compare satellite-based WY with the in-situ-derived WY, and (b) to characterise WY changes in relation to LULC within watersheds according to their corresponding magnitude of LULC. The results of this study will play a crucial role in achieving the United Nations Sustainable Development Goal 6 targets 6.1 and 6.2 through satellite-based analysis of changes in WY due to land use.

2. Materials and Methods

2.1. Study Area

The study was carried out in the Johor River Basin (JRB), located in Johor, the southernmost state of Peninsular Malaysia (PM). The state recorded a total population of about 1,638,219 people at the 2020 census (Department of Statistics Malaysia). The basin covers an area of approximately 2636 km², with a main river length of 123 km [4]. The elevation ranges between 3 and 977 m above mean sea level. The site is an agricultural basin that lies between latitudes 1°30’ and 2°10’N and longitudes 103°20’ and 104°10’E (Figure 1). The primary soil type within the basin is ultisol (Rengam–Jerangau series). This soil is characterised by yellowish–brown sandy clay with moderate permeability and is well-drained, making it suitable for oil palm and rubber plantations [35]. The main land uses of the basin are oil palm and forest [36]. The river flows from the southwest of Johor from south Gunung Belumut at 1010 m elevation. The main tributaries are the Linggiu, Sayong, Penggeli, Jengeli, and Belitong rivers. These tributaries serve as sources of fresh water for the populations of Johor and Singapore. Since the middle of the 1960s, the Public Utility Board (PUB) and the Johor Water Company, Johor, Malaysia, (SAJ) have each drawn approximately 0.25 × 106 km³/day from the JRB [37].

The basin receives an average annual rainfall of 2500 mm, and the rainfall pattern is influenced by two seasons, namely the northeast monsoon (November–February) and the southwest monsoon (May–August). In between these two monsoons, JRB experiences intermonsoon periods that usually happen in March and April and September and October. During the northeast monsoon, most of the eastern coast of Peninsular Malaysia, including JRB, receives heavy rainfall compared to the southwest monsoon and other seasons [38]. The east coast regions receive higher rainfall (>350 mm per month), mainly during December and January.

2.2. Data Used

Five types of data were used in this study. Multitemporal satellite images of Landsat-7-enhanced thematic mapper plus (ETM+) and Landsat operational land imager (OLI) with 30 m spatial resolution, acquired from the United States Geological Survey (USGS) (http://glovis.usgs.gov, accessed on 3 December 2018), were used to generate LULC maps.

Table 1. The attributes of Landsat ETM+ and OLI images were used in the study.

Sensor	Scene ID	Path/Row	Date of Acquisition	* Monsoon
Landsat 7 ETM+	LEO71250582000104EDC00	125/058	29 September 2000	Post-SW
	LEO71250592005356EDC00	125/059	25 September 2005	Post-SW
	LEO71260582010292EDC00	126/058	8 August 2009	SW
Landsat 8 OLI	LE081260592015181EDC00	126/059	2 November 2015	NE

* Notes: SW = southwest, NE = northeast.

summarises the attributes of the four Landsat images used in this study, which were selected based on the longevity of data archiving. An advanced space-borne thermal emission and reflection radiometer digital elevation model (ASTER GDEM), with a spatial resolution of 30 m, was used to derive the vertical slope. A topographic map on a scale of 1:25,000 obtained from the Department of Surveying and Mapping Malaysia (JUPEM) was used to collect ground control points (GCP) for image geometric corrections and the validation of vertical information. TMPA 3B42 version 7 daily data in hierarchical data format (HDF) with a high spatial and temporal resolution from 2000 to 2015 were used to obtain precipitation information for WY extraction. These TMPA data also provide information on rainfall relative error and gauge relative weighting. Because rainfall is the most vital input for modelling WY [39,40], corresponding rainfall observations from rainfall gauges within the JRB were also collected and used for calibration and validation of the TMPA-derived precipitation. A moderate-resolution imaging spectroradiometer (MODIS16A2) satellite data from between January 2000 and December 2015 in HDF, from TERRA and AQUA satellite platforms with 1 km spatial and monthly temporal resolutions, was used to obtain ET. The MODIS is used to provide high-quality ET data [41], containing information on potential latent energy and an ET quality check.

2.3. Methodology

A series of activities were involved in the processing and analysis of all the data acquired for the analysis of changes in JRB. These include LULC mapping, satellite-based water balance extraction, calibration, and validation of TMPA rainfall and MODIS ET data, the extraction of WY changes from individual LULC classes, and analysis of the WY trend. The complete methodological flow chart of this study is illustrated in Figure 2.

2.3.1. LULC Mapping

The LULC mapping was accomplished by Landsat multitemporal mapping to identify and classify different types of land use, such as forests, oil palm, agriculture other than oil palm, urban areas, bare land, and water bodies. This process consists mainly of satellite data preprocessing and image classification.

Pre-Processing of Satellite Images

Usually, deficiencies and errors are found in the raw data obtained via remote sensing satellite-sensor platforms. Therefore, the raw data acquired for the LULC mapping were enhanced by subjecting them to various preprocessing steps to certify their originality. The preprocessing consisted of three major steps: radiometric correction to eliminate systematic errors introduced by sensors and atmospheric conditions, geometric correction to correct for variations in satellite altitude and orientation during data acquisition, and image enhancement to improve the visual quality and highlight important features. All satellite data processing was carried out using the ENVI Digital Image Processing System and ArcGIS system software, as explained in detail below. This differs from image to image depending on the category of information being extracted, the composition of the image scene, and the initial condition of the image.

Initially, the images were radiometrically corrected by applying the rescaling factors and parameters obtained in the metadata file that came with the images. The correction was carried out in two-step processes: (1) the image digital numbers (DNs) were converted to top atmospheric reflectance, and (2) the top atmospheric correction of Landsat OLI and Landsat ETM+ spectral bands was achieved using the FLASH programme of ENVI v.5 software. The atmospheric correction using FLAASH was adopted to lessen the atmospheric uncertainties of nadir-viewing images via inclusivity to correct the effects of adjacency. This is vital to minimise scattering effects [42].

The image-to-map geometric correction technique uses 30 GCP extracted from the topographic maps to correct the images’ [43]. Similarly, corrections were carried out in two-step procedures: (1) a second-degree polynomial was employed to transform the images to map geometry, and (2) the pixel intensity values were created into the transformed geometry by applying the nearest neighbour resampling scheme. The root mean square error (RMSE) of the transformation was ensured to be ±0.5 pixels, and the chosen nearest neighbor resampling scheme in each case was used to avoid loss of details.

Satellite Images Classification

The maximum likelihood classifier, a popular supervised classification method, was used in the LULC mapping. The classification accuracies were evaluated using cross-validation statistics through confusion matrices using a set of 600 stratified random points. An overall accuracy that certifies the minimum threshold of 85% is needed for the effective and steadfast analysis and modelling of the LULC changes [44]. The entire classification assessment comprises the confusion matrices, producers, users, and overall accuracies with the Kappa indexes of the classification images.

2.3.2. Satellite-Based Water Balance Extraction

Satellite-based extraction of water balance operation is focused on the derivation of WY from the rainfall and ET data from TMPA and MODIS satellites, respectively. The WY is normally determined according to the concept of the water balance equation method, which extracts the WY from the equilibrium of precipitation minus ET (see Equation (1));

$$W{Y}_R=P_i-{ET}_{i,}\pm ∆\mathrm{S}$$

(1)

where WY_R is water yield (mm), P_i is precipitation (mm) for the ith month, ET_i is evapotranspiration for the ith month [43], and ∆S is the change in soil moisture, which is insignificant in this study environment [18].

This process was realised by applying digital image processing to obtain the spatial aspect where the WY for each pixel is computed. The WY equation is a standard operation method applied by the Department of Irrigation and Drainage Malaysia (DID). The parameters used in the process of computing the WY are calibrated rainfall and calibrated ET from MODIS. In the computation of WY, there are three components inside the water balance, and all three components cover the initial stages of the water cycle up to the WY. The two main elements that were determined were ET and rainfall.

The rainfall intercepts the LULC, which infiltrates the soil, producing soil moisture; it is in this coincidental phenomenon that ET occurs. Since ET occurs when vegetation exists because of stomatal activities, soil moisture was not addressed. The difference between rainfall and ET was calculated. The results of the two elements can be positive or negative. A positive value shows that precipitation is higher than the ET. This shows the possibility that the WY can be harvested. If the outcome is negative, it means that ET is higher than rainfall.

Before entering rain and ET from TMPA and MODIS satellites into Equation (1), all these data sets were calibrated with the corresponding ground-based observations. Once calibrated, an independent assessment of the calibrated datasets was also performed as part of the validation process. Satellite-based WY assessment was carried out by (1) validation against the observed river flow, (2) comparison with similar studies in nearby watersheds using a water balance equation, and (3) modelling WY using the soil water assessment tool (SWAT).

Modeling WY Using Soil Water Assessment Tool (SWAT)

The SWAT is a semi-dispersed, comprehensive, process-based catchment, and time-uninterrupted model applied to model the potential influence of LULC changes and tue supervision of water quality and quantity [45,46,47,48,49]. SWAT was established by The United States Agricultural Development (USAID) developed the model using its research services in the 1998. The model, combined with a geographical information system (GIS), permits digital LULC, topography, and soil data input. In addition, the model permits the modelling of evaporation and water yield losses of drainage channels. The model repeatedly enhanced SWAT to the very recent version of 2012 [50].

Arnold et al. [49] highlight the model’s flexibility in integrating upland and channel procedures and the modelling of land management. SWAT has several progressive and global applications [50]. The literature has connected to SWAT in several ways [17]. It is an appropriate model for carrying out influence studies, as the effect of LULC, climate change, or both on hydrology has been assessed by applying the SWAT model [51,52]. Additionally, the model has proved its ability to model water changes in regions with inadequate data readiness [53,54]. SWAT is gradually being applied on a larger scale [55,56].

In this study, the model was run in ArcSWAT 2012 interface, which works as an extension in ArcGIS 10.3 software. The SWAT model is freely available on the internet and downloadable from the official SWAT website https://swat.edu/software/arcswat/ (accessed on 1 November 2018). Generally, SWAT modelling is made up of six main steps, as shown in Figure 3. The first step needed during the model run set-up is watershed delineation using a digital elevation model, and the topographic characteristics of the watershed are estimated. The next step is the analysis of hydrological response units (HRU), in which layers of LULC and soil maps were added, and slope ranges were defined. Subsequently, climate station data were added. Next, the input parameters were edited, checked, and validated through data-quality checking, data-sorting or data-processing. After running the model successfully, sensitivity, calibration, and validation were carried out. Similar studies were also carried out [57]. Finally, the extra sub-basin outlet was added based on the location of the hydro-gauging station and

Table 2. Climate parameters in SWAT model database.

Parameters	Details
WLATITUDE	Latitude of the climate station
WLONGITUDE	Longitude of the climate station
WELEV	Height of the climate station in meters above mean sea level
RAIN_YRS	Number of years of maximum monthly half-hour rainfall data used to define values for average per month of the year
TMPMX	Average daily maximum air temperature for each month in degree Celsius (°C)
TMPMX	Average daily maximum air temperature for each month in degree Celsius (°C)
TMPSTDMN	The standard deviation for daily minimum air temperature for each month in degree Celsius (°C)
PCPMM	The average total of monthly precipitation in milometers (mm)
PCPSTD	The standard deviation for daily precipitation for each month expressed as in mm of water per day.
PCPSKW	The skew coefficient for daily rainfall for each month
PR_W1	Possibility of a wet day following a dry day for each month
PR_W2	Possibility of a wet day following a wet day for each month
PCPD	mean numbers of days of rainfall for each month
RAINHHMX	Most extreme half-hour precipitation for each month
SOLARAV	Average daily solar radiation for each month.
DEWPT	Average daily few points temperature per month of the year in degree Celsius (°C)
WNDAV	Average daily wind speed for each month (m/s)

Source: [58].

. The primary purpose of adding an outlet point at the Rantau Panjang Monitoring Station is to compare simulated and observed discharge.

LULC and soil map are inputs to determine the land/soil categories required to establish hydrological response units (HRU). The SWAT categorised the closes climate station to the centroid of sub-catchment and applied it to all the HRUs inside the sub-basin.

Modification of the SWAT model involves the climate and soil parameters in the database. For the climate, SWAT requires daily variables to populate the weather matrix with averages for each month of the year over the total period covered by the station. The parameters are generated from a daily record of the observed data (ideally more than 30 years). For this study, the climate data were from 1980 to 2009, while discharge data were from 1970 to 2018. Table 2 demonstrates the climate parameters that are essential to modify the weather generator (WGN) database in the SWAT model.

Soil information varies from place to place, so modification of the soil database in the SWAT is significant to ensure the better modelling of water within the layers of soil. The soil parameters required to modify the SWAT model’s soil database are presented in

Table 3. Parameters for sensitivity analysis in SWAT model.

Parameter	Details
TEXTURE	Texture of soil layer
HYDGRP	Soil hydrologic group
SOL_ZMX	Maximum rooting depth for soil profile (mm)
ANION_EXCL	A fraction of porosity (void space) from which union are excluded
SOL_CRK	Crack volume potential of soil
SOL_PH1	A soil PH of the first layer of soil.
SOL_Z1	Depth from the soil surface to the bottom of the first layer of soil (mm).
SOL_BD1	Moist bulk density of the first layer of soil (g/cm3)
SOL_AWC1	Accessible water capacity of the soil layer number one (mm)
SOL_K1	Wet hydraulic conductivity of the first layer of soil (mm/h).
SOL_CBN1	The organic carbon content of layer one of the soil (%)
CLAY1	The clay content of layer one of the soil (%).
SILT1	Silt content of layer one of the soil (%).
SAND1	Sand content of layer one of the soil (%)
ROCK1	Rock fragment content of layer one of the soil (%)
SOL_ALB1	Moist soil albedo of layer one of the soil.
NLAYERS	Number of layers in the soil
SOL_EC1	Soil electrical conductivity of the first layer of soil(ds/m).
SOL_CAL1	Calcium carbonate content of layer one of the soil (%)
USLE_k1	USLE equation of soil erodibility (K) factor of the first layer of soil.

Source: [58].

. According to the water balance Equation (5), SWAT carries out the simulation.

$$({SW}_t={SW}_o+\sum _{i=1}^t(R_{day}-Q_{surf}-E_a-W_{seep}-Q_{gw})$$

(2)

where SW_t = final soil water content (mm), SW_o = initial soil water content on day i (mm), t = time (days), R_day = amount of precipitation on day i (mm), Q_surf = amount of surface runoff on day i (mm), E = amount of evapotranspiration on day i (mmH₂ O), W_seep = amount of water entering the vadose zone from the soil profile on day i (mm), and Q_gw = amount of return flow on day i (mm).

The parameters listed in Table 3 are considered the most sensitive parameters that affect the hydrological responses of the river basin. These parameters are applied to control the amount of erosion from the channel and its catchment as they affect the rate of runoff, sediment, and soil nutrient loss to the maximum extent.

The SWAT-Cup 2012 program developed by [59] was applied for a sensitivity analysis of the calibration and validation of the SWAT model. The LULC, DEM, soil and climate data (precipitation, temperature minimum and maximum, solar radiation, relative humidity, and wind speed) were applied as SWAT inputs. The observed monthly discharge (1970–2018) at Rantau Panjang point was applied for calibration and validation. The model was calibrated for 1985–1999 (15 years) after an initial 5-year model warmup period (1970–1984) was used to obtain a better parameterisation of the simulation based on local conditions [49].

The global sensitivity analysis method was used to test 10 parameters with 500 runs (each run has various combinations) performed in parallel with calibration. The fresh parameters obtained during calibration of the model were used for validation of the model. The Sequential Uncertainty Fitting algorithm (SUFI-2), a semi-automatic inverse modelling procedure in the SWAT-CUP, was selected because of its handle and ability to analyse many parameters using the smallest number of model runs [60,61]. A detailed description and processing procedure for the SUFI-2 algorithm within the SWAT-CUP can be found in [62].

2.3.3. Calibration and Validation

Calibration of TMPA Rainfall with Reference to Observed Rain Gauge Data

Monthly rainfall data for nine (9) stations in nine catchments in the JRB from 2000 to 2010 were collected from the Malaysian Meteorological Department (MMD) and DID. Because the reliability of WY measurement depends on the quality of applied rainfall data [63], the performance of the observed rain gauge and TMPA rainfall in the study area was evaluated first using conventional statistical indices [64], using long-term daily rainfall records available at nine catchments in the study area.

The TMPA rainfall observation was also subjected to calibration. This was achieved by a direct comparison of the TMPA observations with the corresponding rainfall obtained from the rain gauge observations using a linear regression analysis approach.

Thus, the calibration function of the annual average TMPA observation is formed and obtained via Equation (3) below:

$${RF}_{calibrated}=C_{RF}+m_{RF}\ast \left({RF}_{TRMM}\right)$$

(3)

where RF_calibrated = calibrated TMPA rainfall, m_RF = slope, c_RF = constant, and RF_TMPA = observed TMPA. The rainfall averages for individual months over 11 years (11 × 12 = 132 months), m_RF and c_RF are 0.93 and 0.04, respectively, with r = 0.90 and R² = 0.73, at s 99% level of significance. These demonstrated that the average monthly data enhance the relationship between the TMPA and the rain gauge data; thus, the average monthly data reduce the differences between the two measurements. The slope and constant are obtained through Equation (3). Various months for each station are suitable for calibrating the TMPA for rainfall at the station level.

Validation of the calibrated TMPA data was also performed. A total of 70% of the rain gauge data randomly selected from nine (9) stations was applied for the TMPA calibration, and the remaining 30% was applied for validation. The performance of the calibrated TMPA was evaluated using RMSE (Equation (4)), where $R_{sat}$ is the calibrated TMPA, $R_g$ is the observed ground rainfall, and n is the number of pixels.

$$RMSE=\sqrt{\left[\frac{1}{n}{\sum }_{i=1}^n{(R_{sat}-R_g)}^2\right]}$$

(4)

Calibration and Validation of MODIS ET Data

The ET obtained from the MODIS satellite data products is the MODIS high-level data product formatted by tiles, and the coordinates of the study area were applied to download the ET. This study used ET retrieved from the digital number (DN) of the MODIS 16A data product. The ET values were multiplied by a constant to convert DN into millimetres per month (Equations (5) and (6)).

$$E_T={MODIS 16A}_{HDF xa}$$

(5)

$${MODIS 16A}_{HDF}=\frac{S x A x C_p\left(\frac{e_{sat}-e}{r_a}\right)}{s+y(l+\frac{r_s}{r_a})}$$

(6)

where ET is the total ET estimates at monthly intervals (mm/month), MODIS16A_HDF is unitless ET in HDF format, a is a constant which is set to 0.1, s = d($e_{sat}$)/dts = denotes the slope of the curve relating to saturated water-vapour pressure ($e_{sat}$) to temperature, A is the available energy divided between sensible heat, latent heat, and soil flux on the land surface, ρ is the air density, $C_{\mathsf{\rho }}$ is the specific heat capacity of air, $r_{\mathrm{a}}$ is the aerodynamic resistance, and $r_{\mathrm{s}}$ is the surface resistance. Surface resistance was parameterized by applying the satellite leaf area index and vegetation fraction cover. γ is the psychrometric constant. The detailed procedure for estimating MODIS 16A (ET) is presented in Figure 4.

Satellite-based ET provides data on a global scale. Although there are numerous validations of those data sets for various climates, including in America and Asia, specific validation, especially for the tropical region, was carried out by [18]. These studies suggested that it is necessary to calibrate the data. In addition, errors related to seasonal variation have been shown to require calibration by using a linear regression function (Equation (7)) based on monsoon characteristics.

$${ET}_{calib}=a\ast {ET}_{MODIS}+c,$$

(7)

where ET_calib is the calibrated satellite ET and a and c are the calibration coefficients with values of 0.36 and 54.7, respectively.

The calibration of the ET on MODIS 16A2 was performed by computing the ET from the Kluang meteorological station [18], which is the nearest weather station, using the Penman–Monteith (PM) method [65], involving surface and water-vapor aerodynamics based on Equation (8):

$$ET=\frac{0.408\Delta \left(R_n-G\right)+\gamma \frac{900}{T+273}{u}_2(e_s-e_a)}{\Delta +\gamma (1+0.34u_2)}$$

(8)

where $ET$ = reference ET rate (mm month⁻¹), $\mathsf{\Delta }$ = slope of the vapor pressure curve (kpa°C⁻¹), $\gamma$ = psychrometric constant (kpa ${℃}^{-1}$) $R_n$ = net radiation at the crop surface (MJ m⁻² d⁻¹), G = soil heat flux density (MJ m⁻² month⁻¹) (scale G is daily assumed to be zero), $e_s$ = saturation vapor pressure (kpa), $e_a$ = actual vapor pressure (kpa), $e_s-e_a$ = saturation vapor pressure deficit (kpa).

The accuracy of the calibrated ET MODIS 16A was checked and validation tests were also conducted. The MODIS data were calibrated using rain gauge data collected in the periods between 2000 and 2006 and validated using data collected in the periods between 2007 and 2010. The performance of the calibrated MODIS 16A was evaluated using RMSE. The results are shown in Section 3, Results.

2.3.4. Satellite-Based Water Yield Changes from Individual LULC Classes

The extraction of changes in total WY from individual LULC classes was carried out using a logistic function to establish the correlation between the WY from the water balance equation and those derived from satellite images of different types of LULCC. The developed models used in estimating WY from the LULC classes are detailed siEquation (9):

$${LULC}_{WY}=AOOP+BL+FRST+OP+UB+WB,$$

(9)

where:

• ${LULCY}_{TWY}$ = total WY from LULC classes;
• AOOP = WY from agriculture LULC class;
• BL = WY from bare land class;
• FRST = WY from forest class;
• OP = WY from oil palm class;
• UB = WY from urban class;
• WB = WY from water body class.

Therefore, the total WY of the study area was extracted from Landsat ETM+ and OLI using the developed model. This derivation process was carried out by extracting satellite-based WY (per pixel) from six LULC classes (agriculture other than oil palm, bare land, forest, oil palm, urban and water bodies). Therefore, to understand the developed total WY equation, the coefficients of the variables (which represent the individual LULC classes) are used as an indicator of the contribution of each LULC to the total WY.

The effects of LULC changes on WY are analysed at 5- and 10-year intervals, focusing on LULC changes in 2000–2005, 2005–2010, 2010–2015, and 2000–2015. The LULC changes to total WY ($∆WY)$ for the entire JRB from the changes in LULC classes were analysed using multiple regression analysis, given by (Equations (10)–(13)):

$$\begin{gathered} {∆WY}_{2000-2005}=b_1∆{WY}_{AOOP}+b_2∆{WY}_{BL}+b_3∆{WY}_{FRST}+b_4∆{WY}_{OP} \\ +b_5∆{WY}_{UP}+b_6∆{WY}_{WP} \end{gathered}$$

(10)

$$\begin{gathered} {∆WY}_{2005-2010}=b_1∆{WY}_{AOOP}+b_2∆{WY}_{BL}+b_3∆{WY}_{FRST}+b_4∆{WY}_{OP} \\ +b_5∆{WY}_{UP}+b_6∆{WY}_{WP} \end{gathered}$$

(11)

$$\begin{gathered} {∆WY}_{2010-2015}=b_1∆{WY}_{AOOP}+b_2∆{WY}_{BL}+b_3∆{WY}_{FRST}+b_4∆{WY}_{OP} \\ +b_5∆{WY}_{UP}+b_6∆{WY}_{WP} \end{gathered}$$

(12)

$$\begin{gathered} {∆WY}_{2000-2015}=b_1∆{WY}_{AOOP}+b_2∆{WY}_{BL}+b_3∆{WY}_{FRST}+b_4∆{WY}_{OP} \\ +b_5∆{WY}_{UP}+b_6∆{WY}_{WP} \end{gathered}$$

(13)

where:

• WY_2000–2005 = changes total WY between 2000 and 2005;
• b₁, b₂, b₃, b₄, b₅, and b₆ are coefficients of WY changes due to the LULC classes of agriculture, bare land, forest, oil palm, urban, and water body, respectively.

The assessment of changes in WY for the changes in LULC changes was carried out, focusing first on the model-fitness (based on coefficients and significance level) on WY from each LULC class, and secondly on the validation (using RMSE) of these LULC WY.

2.3.5. Mann–Kendall Statistics for Trends Analysis

Mann–Kendall (MK) statistics [63] were employed to analyse the trends of the WY (decreasing or increasing) within the JRB. The MK test was selected because it has the following advantages: (i) it does not require normally distributed data, (ii) it has been widely used by the World Meteorological Organisation (WMO) and (iii) it can handle outliers. In the trend test of Mann–Kendall, individual data in a series are compared to all succeeding data in the series. In this method, the changes between each consecutive value are computed to show increasing (+1), decreasing (−1), and neutral (0) signs. The MK (S) for a given data series x₁, x₂, x₃….., and xn (LULC WY) were computed using Equations (14)–(17):

$$S={\sum }_{k=1}^{n-1}{\sum }_{i=k+1}^{n}sign(x_1-x_k)$$

(14)

where

$$sign\left(x_1-x_k\right)=\left\{\begin{array}{c}+1if\left(x_1-x_k\right)>0\\ 0if\left(x_1-x_k\right)=0\\ -1if\left(x_1-x_k\right)<0\end{array}\right.$$

(15)

where n is the number of data in the series while $x_1$ and $x_k$ are the ranks for the data. The significance of the trends in the data series can be calculated using statistics Z in Equation (15), as below:

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}}ifS>0\\ 0\hspace{1em}ifS=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}}ifS<0\end{array}\right.$$

(16)

where Var (S) is the variance in S positive and negative Z values, indicating the direction of the trend that exists in the time series. In this expression, V(S) variance, and in case of tied ranks it is given as in Equation (16):

$$V(S)=\frac{n(n-1)(2n+5){\sum }_{i=1}^m{t}_1(t_1-1)(2t_1+5)}{18}$$

(17)

where m is the number of tied groups and ${\mathrm{t}}_1$ is several observations in the ith group. If there is no dependence, then V(S) can be obtained through Equation (16).

3. Results

There are two main results in this study. These include a comparison of satellite-based WY with in-situ-derived WY and a characterization of WY changes in relation to LULC within watersheds according to their corresponding LULC magnitude. However, validation and calibration results are first highlighted as valuable additions to further enhance the findings.

3.1. Calibration and Validation of Satellite Data

3.1.1. Calibration and Validation of TMPA

Table 4. Evaluation of the performance of the monthly rainfall of TMPA versus the observation of the rain gauge in annual rainfall estimation at all nine catchments in the study area.

S/NO.	Stations	Lat. (N)	Long. (E)	Rain Gauge Obs.	TMPA Obs.
				Annual Rainfall (mm/yr)	R2	NSE	BIAS%
1	Sg. Jengeli	01°57′00″	103°39′00″	2268.38	0.67	0.65	9
2	Sg. Johor	01°45′30″	103°50′00″	2158.92	0.65	0.63	10
3	Sg. Johor *	01°35′30″	103°56′30″	2592.13	0.69	0.67	8
4	Sg. Layang	01°32′30″	103°53′00″	2106.36	0.64	0.62	−9
5	Sg. Lebak	01°49′00″	103°48′00″	2323.95	0.73	0.71	−7
6	Sg. Linggiu	01°59′30″	103°40′30″	2085.75	0.60	0.60	−8
7	Sg. Sayong	01°52′30″	103°30′00″	2435.79	0.75	0.73	8
8	Sg. Seluyut	01°45′00″	104°00′00″	2388.05	0.70	0.69	4
9	Sg. Semangar	01°44′00″	103°40′00″	2487.75	0.68	0.66	9

* Upper Sg. Johor.

presents the monthly observed rainfall data from nine rain-gauge locations compared with TMPA rainfall data. Three standard statistical indices were obtained, namely: coefficient of determination (R²), Nash–Sutcliff efficiency (NSE) and Bias were obtained to evaluate the performance of the observed rainfall and TMPA data. They showed good agreement with TMPA monthly rainfall at the nine observed locations in JRB, with R² in the range of 0.60–0.75, NSE values > 0.65 at six out of nine locations, and bias at less than 10% at all the locations.

Figure 5 summarises the linear regression analysis used for the calibration of TMPA, showing the relation between the observed TMPA and the corresponding rainfall gauge data from 2000 to 2010, from the monthly rainfall of the nine (9) stations, obtained for 132 months in 11 years. The solid line displays the linear relationship, while the dotted line is a control indicator y = x.

Figure 6 presents a comparison of uncalibrated and calibrated TMPA data using RMSE. The calibration minimises the random errors inherent in the TMPA observation, reducing the RMSE by up to 16.3%. The RMSE rainfall values estimated from TMPA calibrated data indicate that the rainfall data sets improved after calibration in terms of the RMSE values. This data calibration improvement pattern agrees with the temporal variation, where the difference between calibrated TMPA and ground measurement serves as evidence.

3.1.2. Calibration and Validation of MODIS ET

The calibration function derived from the regression of the ET MODIS 16A product against the corresponding ground ET at the Kluang station is presented in Figure 6. This study achieved accuracy using RMSE (±43 mm) and the determination coefficient (R² = 0.642). These assessments were realised using eleven (11) years’ monthly mean of both satellite and in situ observations (Figure 7). RMSE plots for MODIS 16A ET and data from the Kluang meteorological station in 2007 and 2010 in various seasons, showing changes before and after the calibration, are presented in Figure 8. The variation in error ranged from 16 to 26 mm in the uncalibrated and calibrated MODIS ET. The calibration of MODIS ET led to improvements in WY monitoring, similarly to the research conducted by [18].

3.2. Assessment of Fully Satellite-Based Water Yield Compared to In-Situ-Derived Water Yield

3.2.1. Validation with River Flow Data

The satellite-based WY for JRB was also analysed to compare flood and drought occurrences. The time series graph of WY versus flood and drought is shown in Figure 9, where this catchment was reported to experience 100-year flooding in recent years [66]. The peak and lowest WY derived from this satellite-based method correspond well to the records of local flood and drought occurrences in JRB, respectively. The satellite-based and in-situ-derived WY is shown in

Table 5. The satellite-based and in-situ-derived water yield.

S/No.	Months	WY from Satellite-Based Water Balance Equation (mm/Month)				WY from SWAT (mm/Month)
		2000	2005	2010	2015	2000	2005	2010	2015
1	January	165.76	175.52	124.29	215.86	169.23	135.93	157.87	200.67
2	February	119.49	233.95	250.14	128.24	103.88	247.23	226.17	122.59
3	March	103.86	100.66	118.17	161.44	139.47	55.66	126.02	173.01
4	April	122.57	85.26	94.83	133.98	143.26	42.46	100.00	117.4
5	May	82.38	79.35	62.40	97.51	88.06	83.47	69.59	160.63
6	June	35.17	43.82	54.49	70.20	38.08	34.09	72.12	82.49
7	July	64.88	64.98	66.97	96.81	76.48	82.08	103.29	110.53
8	August	63.60	103.94	75.91	86.81	83.36	118.55	88.80	93.97
9	September	77.49	82.92	84.60	122.44	72.36	124.4	77.35	146.65
10	October	65.63	171.61	86.52	254.29	78.96	189.34	79.73	229.29
11	November	94.40	133.05	294.25	230.85	83.35	174.51	256.78	217.80
12	December	228.87	161.15	219.50	367.89	194.90	171.16	155.32	317.74
Total WY (mm/yr)		1224.10	1436.21	1532.07	1966.32	1251.41	1458.88	1513.04	1972.77

.

3.2.2. Validation with Soil Water Assessment Tool (SWAT)

The observed and modelled monthly discharge at the Rantau Panjang discharge point is shown in Figure 8. The calibration period started from January 1985 to December 1999, while the validation period was from January 2000 to December 2009. The NSE, R² and PB values were (0.67, 0.67) and (−3.1), respectively, for the calibration period, whereas (0.63), (0.65) and (−1.9) are the values of NSE, R² and PB, respectively, for the validation period. Based on [67], for the calibration period, the NSE values demonstrate that the SWAT model for the JRB was deemed to have too good a performance for the calibration and validation times.

In Reference [68], the model demonstrates good discharge modelling performance in the Bukit Merah Reservoir, Malaysia, with an R² of 0.87 and 0.69 for the calibration and validation periods, respectively. For the NSE method, the SWAT output is 0.79 for calibration and 0.60 for validation periods. The performance of the SWAT model was considered ‘very good’ in the calibration period and ‘satisfactory’ in the validation period. Figure 10 demonstrated the observed and simulated means of monthly discharge recorded at the Rantau Panjang station.

The performance of the SWAT model for the period of calibration is better than that of validation; this might be because the temporal variations in the SWAT model parameters were not effectively considered. In addition, the hydrograph validation period demonstrates an overestimation of the simulated discharge during the southwest monsoon season, which was also found in other studies [66]. This may be because of the occurrence of extreme flood, where the model poorly matched the peak flow. The application of the model in Malaysia is still limited; hence, SWAT model calibration and validation in this study demonstrate that it is a consistent tool for hydrology cycle modelling in Peninsular Malaysia.

3.2.3. Comparison of Similar Studies

The WY estimates for selected catchments in Malaysia, which have quite similar characteristics in terms of land use and rainfall, are compared with the present study (

Table 6. Estimate of water yield values from selected basins in Malaysia.

Watershed/Catchment	Type of Watershed/Catchment	Size (km2)	Mean Annual Precipitation (mm)	Total Water Yield (mm yr−1)	References/Approach Satellite
Johor River	Oil palm plantation, forest, and semi-urban	2367.17	2500	710	This study/Satellite
Langat River	Semi-Urban	1257.70	2401	1207	[18]/Satellite
Layang	Semi-Urban	33.61	2690	1334	[69]/SWAT/Satellite
Pendang Terap	Forest	1032.30	2406	868	[18]/Satellite
Hulu Perak	Forest	857.30	2641	687	[18]/Satellite
Kenyir Lake	Forest	1260.00	2606	1473	[18]/Satellite
Johor River	Semi-Urban	2636.50		788	[18]/Satellite
Johor River	Semi-Urban	1652.00	2500	706	[29]/SWAT/Satellite
Hulu Langat	Forest	390.26	2453	742	[70]/SWAT/Satellite

). Most of the analyses used the satellite to estimate the WY, except for the studies by [29] for the Johor River and by [69] for the Layang River, which combined SWAT and satellite data. The annual WY for various catchments ranges from 706 mm/year for JRB [29] to 1473 mm/year for Kenyir Lake [18]. For JRB, the WY estimates found in the earlier studies are quite similar to those in the present study.

3.3. Analysis of LULC Changes and Water Yield

Figure 11 presents the analysis of the change in LULC and WY in JRB between 2000 and 2015. Figure 11a depicts the spatio-temporal pattern of LULC classes, while Figure 11b presents their corresponding WY. The average WY for the entire catchment is presented in Figure 11c. The LULC distribution (

Table 7. Individual LULC classes’ area water yield estimated for 2000, 2005, 2010, and 2015, as well as JRB WY for yearly comparison.

LULC Classes	Yr 2000		Yr 2005		Yr 2010		Yr 2015
	Area	WY	Area	WY	Area	WY	Area	WY
	(ha)	(mm/yr)	(ha)	(mm/yr)	(ha)	(mm/yr)	(ha)	(mm/mr)
Agriculture other than oil palm	11,943.5	535.28	24,412	717.18	33,480.8	879.59	38,169.3	983.06
Bare land	2885.4	856.7	2409.66	632.46	2145.87	314.88	5246.1	306.04
Forest	101,717	77.33	78,060.24	196.29	46,648.6	277.53	25,624.3	358.76
Oil palm	95,051.1	101.03	103,527	118.92	108,755.6	237.97	123,164	280.4
Urban	10,683	170.95	15,517.1	416.41	33,091	549.4	33,924.1	592.9
Water bodies	14,437.4	63.37	12,791.4	55.7	12,595.5	53.08	10,589.6	29.66
Total Water Yield (mm/yr) In JRB		606.32		756.21		803.07		1074.1

) shows that oil palm and agriculture, other than oil palm, were the main land-cover classes, followed by urban, while there was a decline in the forested area. The accuracy of the LULC classification, summarised in

Table 8. Classification accuracy of LULC maps of Johor River Basin.

LULU Categories	2000 (%)		2005 (%)		2010 (%)		2015 (%)
	Prod.	User	Prod.	User	Prod.	User	Prod.	User
Agricultural area	86.33	92.64	91.33	85.50	86.33	87.83	92.00	89.46
Bare Land	89.67	89.67	88.00	91.10	91.33	92.92	90.26	87.55
Forest	93.00	88.48	85.00	85.93	98.00	88.91	90.68	92.45
Oil Palm	94.67	85.88	88.00	86.52	86.33	86.33	88.76	89.77
Urban	88.00	92.74	86.33	87.83	91.33	92.92	89.44	90.26
Waterbody	91.33	94.55	89.67	89.67	88.00	92.74	92.33	89.84
Overall Kappa	89.00		85.67		88.67		88.05
Overall Accuracy	90.50		87.72		90.22		86.89

, revealed high accuracy, which could not be unrelated to the use of only six classes in the classification. However, the results of the accuracy assessment show that the image classifications met the acceptance criteria.

The study found that between 2000 and 2015, WY from agricultural land other than oil palm, forest land, oil palm, urban areas, and water bodies increased by 20.58%, 12.93%, 8.24%, 19.39% and 1.55%, respectively, while WY from bare land decreased by 25.31%. Mann–Kendall statistics show that WY in JRB with respect to LULC showed an increasing trend between 2000 and 2005 but a decreasing trend between 2010 and 2015 (Figure 11b). The sub-basin with the highest average WY was Sg. Layang (160.32 mm/month) while the sub-basin with the lowest was Sg. Sayong (131.95 mm/month) (Figure 11c).

The results of the WY for all the catchments showed that the WY increased in 2015 compared to 2000, particularly for forest (Figure 12) and agriculture (other than oil palm). Oil palm and urban WY also increased, while bare land and water bodies WY decreased in their respective catchments from 2000 to 2015. However, the WY of individual LULC classes, estimated using a water balance equation in JRB for 2000, 2005, 2010, and 2015, respectively, is presented in

Table 9. Individual LULC classes’ water yield estimation using a water balance equation in JRB for 2000, 2005, 2010, and 2015.

S/No.	LULC Classes	2000 WY (mm/yr)	2005 WY (mm/yr)	2010 WY (mm/yr)	2015 WY (mm/yr)	Δ 2000–2015 (%)	p-Val.	RMSE
1	Agriculture other than oil palm	535.28	717.18	879.59	983.06	20.58	0.001	0.64
2	Bare land	856.7	632.46	314.88	306.04	25.31	0.001	0.58
3	Forest	77.33	196.29	277.53	358.76	12.93	0.001	0.56
4	Oil palm	101.03	118.92	237.97	280.4	8.24	0.001	0.6
5	Urban	170.95	416.41	549.4	592.9	19.39	0.001	0.62
6	Water body	63.37	55.7	53.08	29.66	1.55	0.001	0.59

.

3.4. Characterisation of WY in All Watersheds

The variation in WY intensities within LULC is shown in Table 8, and the changes in the LULC class versus WY changes are shown in Figure 13. Consequently, between 2000 and 2015, agricultural land increased by 18.78%, which led to an increase in WY by 20.58% in the same period. The bare land is reduced by 0.19%, leading to a reduction in WY of 25.31%. In contrast, a reduction in forest area to 33.40% increased the WY by 12.93%. Furthermore, an increase in oil palm area to 13.51% increased the WY by 8.24%, whereas an increase in urban areas to 18.50% increased the WY by 19.39%. The water bodies are reduced to −1.72%, while the WY extracted from the water bodies is also reduced by 1.55% (

Table 10. LULC changes 2000–2015 and WY changes.

* LULC Classes	LULC Changes 2000–2005		LULC Changes 2005–2010		LULC Changes 2010–2015		LULC Changes 2000–2015
S/No.	Ha	%	Ha	%	Ha	%	Ha	%
1	12,455.73	5.26	9081.63	3.84	22,838.76	9.67	44,355.87	18.78
2	−482.94	−0.20	−256.59	−0.11	293.76	0.12	−448.11	−0.19
3	−23,590.98	−9.97	−31,476.87	−13.30	−23,920.47	−10.13	−78,897.78	−33.40
4	8429.49	3.56	5261.31	2.22	18,265.68	7.73	4644.45	13.51
5	4824.90	2.04	17,583.48	7.43	21,307.23	9.02	43,695.90	18.50
6	−1636.20	−0.69	−192.96	−0.08	−2253.60	−0.95	−4061.43	−1.72
Changes in Water yield (mm)	149.89	16.02	46.86	5.01	271.03	28.97	467.78	50.00

* Notes. 1 = agriculture other than OP; 2 = bareland; 3 = forest; 4 = oil palm; 5 = urban; and 6 = water body.

).

As outlined in Equation (1), in modelling the total WY attributed to LULC changes, it was found that the use of the AIC was the lowest at −394.78, suggesting the model as the optimum, as in Equation (18), and the full solution tabulated in

Table 11. The equation to best predict water yield.

LULC Classes	Estimate	Std. Error	t Value
Intercept	−0.004165 ***	0.001056	−3.942
Agriculture	1.643754 ***	0.040208	40.881
Bare land	1.093582 ***	0.060279	18.142
Forest	1.003378 ***	0.089682	11.188
Oil Palm	0.688391 ***	0.149676	4.599
Urban	0.752166 ***	0.016001	47.008

*** indicates a significant difference at p < 0.001.

. The AIK is a statistical method to evaluate how well a model fits the data from which it was generated.

However, the water body does not contribute significantly to WY. The strongest parameter that contributed to changes in WY was agriculture. This could be supported by the fact that the water balance equation used is based on the runoff model; hence, it is evident of there being no significant runoff over water bodies.

$$y_i=-0.004165+1.644x_1+1.094x_2+1.033x_3+0.688x_4+0.752x_5+{\epsilon }_i$$

(18)

where yi is total water yield (in mm unit) for watershed; x₁, x₂, x₃, x₄, x5 is amount for WY (unit in mm) for land cover classes of agriculture (other than oil palm), bare land, forest, oil palm and urban, respectively, and $\epsilon is error such that \epsilon ~N(0,0.037)$.

However, there was no viable solution for Equation (17) in a further investigation of which LULC class changes affect the intensity of respective WY changes due to “limited LULC changes.” The LULC changes do not show many changes within a short period, such as a less than a 5-year interval, except for the small percentage of urbanised areas. From 2000 to 2015, only four LULC maps were generated by classifying Landsat TM and OLI data acquired in 2000, 2005, 2010, and 2015. On the other hand, the average monthly satellite-based WY could be obtained from the daily WY derived from satellite precipitation and ET rates, producing adequate dependent parameters.

4. Discussion

Changes in WY have become a crucial issue in sustainable development across the globe. One of the obvious reasons for this change in WY is the LULC change. Satellite-based earth observations, as an important source of data for several aspects of the earth, have the key benefits of a synoptic view of the earth surface, regular and repeatable observation, monitoring of remote and inaccessible areas, and time series observation. Therefore, they can offer valuable data that can be applied to precisely estimate WY [8,71] and its changes over time due to LULC changes.

The analysis of WY in JRB between 2000 and 2015 revealed that changes in land use can have a significant impact on WY, as is the case in some forested ecosystems in the United States of America (USA) [72] due to the impact of landcover changes on evaporation, streamflow, and runoff [72,73,74]. However, in line with one of the recommendations by [72] to consider other landcovers, this study revealed that the highest increase in WY over the study area came from agricultural land other than oil palm, which increased by 18.78% and resulted in a 20.58% increase in WY. An increase of 13.51% in oil palm area also resulted in an 8.24% increase in WY. This demonstrates that an increase in agricultural activity can result in an increase in water resources and is likely due to a reduction in infiltration due to the replacement of natural vegetation cover with agricultural products and an increase in surface runoff. The increase of 12.93% in WY due to the reduction in forest area to 33.40% further emphasises the argument of an increase in WY due to an increase in agriculture, which could be due to the conversion of natural vegetation cover to agricultural fields and a reduction in infiltration. Nevertheless, the case of an increase of 19.39% WY due to an 18.50% increase in urban land cannot be unrelated to the increased compaction of the surface, reduced infiltration, and increased surface runoff, as highlighted in [75]. Finally, the 1.55% decrease in WY due to the reduction in water bodies is not surprising given the increased conversion of natural land for agriculture and the reduction in forest, which can expose shallower surface water to scorching evaporation, and, of course, lead to increased consumption due to increased urbanisation.

However, while the increase in water yield found in this study is a positive result, it is critical to recognise the potential negative consequences of natural vegetation conversion to agriculture. Land cover changes in the form of damage to endangered native vegetation can have a wide range of environmental implications [76], which should be carefully evaluated. The possible loss of biodiversity and ecosystem services offered by natural vegetation is one of the major drawbacks of this land cover change. Forests, for example, are home to many unique plant and animal species that contribute to the general health and stability of ecosystems. There is a risk of habitat fragmentation, loss of species richness, and the disruption of ecological processes when forests are converted to agricultural land. These ecological alterations may have long-term consequences for ecosystem functioning, such as diminished pollination, nutrient cycling, and natural pest control, among other things [77,78]. Furthermore, the conversion of natural vegetation to agriculture frequently involves the use of intensive farming practices, such as the application of artificial fertilisers and pesticides. The improper management of these inputs can lead to soil deterioration, water contamination, and a loss of soil fertility over time. Furthermore, the removal of vegetation cover can enhance soil erosion, as already reported in the study area by [79], resulting in sedimentation in surrounding water bodies and potentially detrimental effects on water quality, as reported in [80]. It is also important to think about the social and economic implications of land cover changes. The conversion of natural vegetation to agricultural land may have consequences for local communities, particularly those that rely on forest resources for a living. Some of the social consequences that should be considered include the displacement of indigenous or marginalised communities, the loss of traditional knowledge, and changes in cultural practices.

Notwithstanding, this study highlights the importance of natural ecosystem conservation and minimising land conversion for other uses, as well as the importance of proper water body management in order to maintain their contributions to water resources. Overall, these findings indicate that careful consideration of LULC changes in the JRB is critical to ensure the basin’s water resource sustainability and contributes to the development of a baseline for the effects of LULC changes on WY in the basin, which is highly significant for efficient land-use planning, ecological restoration, and management and guidance for regional socioeconomic development. Furthermore, the study can contribute significantly to related industries, such as the Public Utility Board, in deciding water resource assessment with regard to the quantity of water. The study contributes substantial support to understanding the degree of expansion in urban planning and the distribution of resources within the watershed.

5. Conclusions

In summary, the study reveals satellite-based techniques’ ability to serve as a powerful tool to monitor the effects of spatiotemporal trend mapping and the estimation of LULC changes on WY. This effort was carried out between 2000 and 2015 in the JRB of Peninsular Malaysia. The effects of LULC changes on WY were achieved using a water balance equation technique that determines WY from the equilibrium of precipitation minus ET. The process was achieved by applying digital image processing to obtain WY from individual LULCs for the entire catchment. The method was validated using in situ measurement and derived values and reported good agreement with ground-based rainfall, ET, and river discharge. Between 2000 and 2015, agricultural areas other than oil palm increased to 11.07%, forest decreased to 32.15%, oil palm increased to 11.88%, urban areas increased to 9.82%, and WY increased to 15.76%. The findings will provide valuable information for water resource management and planning, as well as aiding policymakers in taking proactive steps to mitigate the effects of changes in LULC to WY. Consequently, the outcome will fast-track the attainment of the 2030 agenda of United Nations Sustainable Development Goal 6, targets 6.1 and 6.4.