A Multi-Variable Sentinel-2 Random Forest Machine Learning Model Approach to Predicting Perennial Ryegrass Biomass in Commercial Dairy Farms in Southeast Australia

: One of the most valuable and nutritionally essential agricultural commodities worldwide is milk. The European Union and New Zealand are the second-and third-largest exporting regions of milk products and rely heavily on pasture-based production systems. They are comparable to the Australian systems investigated in this study. With projections of herd decline, increased milk yield must be obtained from a combination of animal genetics and feed efficiencies. Accurate pasture biomass estimation across all seasons will improve feed efficiency and increase the productivity of dairy farms; however, the existing time-consuming and manual methods of pasture measurement limit improvements to utilisation. In this study, Sentinel-2 (S2) band and spectral index (SI) information were coupled with the broad season and management-derived datasets using a Random Forest (RF) machine learning (ML) framework to develop a perennial ryegrass (PRG) biomass prediction model accurate to +/ − 500 kg DM/ha, and that could predict pasture yield above 3000 kg DM/ha. Measurements of PRG biomass were taken from 11 working dairy farms across southeastern Australia over 2019–2021. Of the 68 possible variables investigated, multiple simulations identified 12 S2 bands and 9 SI, management and season as the most important variables, where Short-Wave Infrared (SWIR) bands were the most influential in predicting pasture biomass above 4000 kg DM/ha. Conditional Latin Hypercube Sampling (cLHS) was used to split the dataset into 80% and 20% for model calibration and internal validation in addition to an entirely independent validation dataset. The combined internal model validation showed R 2 = 0.90, LCCC = 0.72, RMSE = 439.49 kg DM/ha, NRMSE = 15.08, and the combined independent validation had R 2 = 0.88, LCCC = 0.68, RMSE = 457.05 kg DM/ha, NRMSE = 19.83. The key ﬁndings of this study indicated that the data obtained from the S2 bands and SI were appropriate for making accurate estimations of PRG biomass. Furthermore, including SWIR bands signiﬁcantly improved the model. Finally, by utilising an RF ML model, a single ‘global’ model can automate PRG biomass prediction with high accuracy across extensive regions of all seasons and types of farm management


Introduction
Dairy farms in Australia are largely grazing-based systems (almost 96%), utilising pasture as a nutrient-dense and cheap feed source for dairy animals [1,2]. This precious resource can be efficiently managed if measured, but in busy farming systems such as dairy, it is difficult to observe the whole farm all the time adequately using traditional on-farm pasture biomass measuring approaches. Fulkerson et al. [3] found that feed An extensive array of studies has demonstrated processing options for multispectral sensor bands, such as single bands [35], commonly reported spectral indices, linear spectral unmixing [26], spectral mixture analysis [36], and using the principles of radiation use efficiency [37], and physically-based radiative transfer modelling [38]. Vegetation indices are the most common approach, with the normalised difference vegetation index (NDVI) being the most widely used and reported over the last 50 years and commonly used as a reference index on which to gauge improvement [26,39,40]. Although NDVI has been proven to be an excellent predictor of biomass in many studies, the challenges of fastgrowing high-volume wet canopies, such as those found in dairy pastures, suggest the investigation of a suite of spectral indices (SI). Numata et al. [36] found better correlations with fresh weight than dry weight using two indices that use SWIR bands. SWIR bands have been found to improve the correlations between the measured and predicted dry weight due to sensitivity to leaf water concentration at high biomass (>3000 kg DM/ha). Biomass above 3000 kg DM/ha is a level that commonly exceeds a Leaf Area Index (LAI) greater than 3. LAI > 3 is known to saturate common vegetation indices such as the Normalised Difference Vegetation Index (NDVI) [41]. Broge and Leblanc [42] demonstrated that the effectiveness of SI depends not just on the bandwidth or combination but also on the target's range and environmental factors such as soil classification and colour variations in the soil in both dry and wet conditions and moisture present in the air that affect that range.
With the ever-increasing data from modern satellite sensors, selecting a proper empirical modelling framework is essential in landscape research. The linear regression analysis approach has been favoured for data analysis, particularly for satellite-based work [31,32,43]. Guerini Filho et al. [32] considered multiple linear regression analyses for natural grassland biomass estimation by SI and obtained a coefficient of determination (R 2 ) less than 0.70. Although the accuracies were satisfactory, the major disadvantage of the approach is the site-dependent research and inability to capture the complex non-linear patterns found in the data. Therefore, the algorithm will likely have greater errors for the diverse pasture types in dairy management systems due to the design of the study presented by Guerini Filho et al. [32].
Machine Learning (ML) algorithms have efficiently handled large datasets, particularly the ones with non-linear patterns. Recent progress in digital agriculture has also triggered the necessity of state-of-the-art ML approaches to increase the efficiency of model prediction with further automation [44]. ML has been used in agricultural applications, particularly for pasture biomass prediction [25]. However, there is no universal measure of ML suitability [45]. Among different ML techniques, non-parametric supervised classifiers such as support vector machine (SVM) [46], artificial neural networks (ANN) [47,48], k-nearest neighbor (kNN) [49], classification and regression tree (CART), and supervised techniques such as Random Forest (RF) have been widely used for predicting biomass based on satellite images [50,51].
The significant objectives of this study are summarised as the following: • Combine S2 and ML modelling approaches to predict perennial ryegrass biomass across a range of regions and seasons with an accuracy of +/−500 kg DM/ha. • Provide technical evidence that utilising SWIR bands can improve the ability to predict pasture yields above 3000 kg DM/ha and, therefore, enable measurements of highyielding pastures at any stage in the growth cycle in irrigated and dryland farm management systems. • Examine the pasture biomass prediction model quality through a fusion of S2 sensorderived datasets and broad management and seasonal datasets. • Show that it is possible to predict pasture biomass across large regions and growing seasons on commercial dairy farms with one 'global' model with an extensive ground sampling campaign and the use of numerous bands and SI of the S2 sensor and the ML modelling framework.

Study Site Location, Soil, and Climate and Sampling Design
A total of 11 commercial dairy farms were part of this study. Seven farms were intensively sampled and contributed to the calibration and validation dataset. The remaining farms were sampled once to obtain fully independent validation data. Four farms were in the Northern Irrigation Region of Victoria, two in the Macalister Irrigation District in far eastern Victoria, two in West Gippsland, two in the southwest of Victoria, and one in the southeast of South Australia, as illustrated in Figure 1. Four of the eleven farms were dryland farms, and the remaining properties applied irrigation water to some or all of the grazed pastures using flood irrigation or overhead sprinklers. Farm sizes ranged from 72 to 970 ha.
Remote Sens. 2023, 15,2915 4 of 3 • Show that it is possible to predict pasture biomass across large regions and growin seasons on commercial dairy farms with one 'global' model with an extensive groun sampling campaign and the use of numerous bands and SI of the S2 sensor and th ML modelling framework.

Study Site Location, Soil, and Climate and Sampling Design
A total of 11 commercial dairy farms were part of this study. Seven farms were inten sively sampled and contributed to the calibration and validation dataset. The remainin farms were sampled once to obtain fully independent validation data. Four farms were i the Northern Irrigation Region of Victoria, two in the Macalister Irrigation District in fa eastern Victoria, two in West Gippsland, two in the southwest of Victoria, and one in th southeast of South Australia, as illustrated in Figure 1. Four of the eleven farms were dry land farms, and the remaining properties applied irrigation water to some or all of th grazed pastures using flood irrigation or overhead sprinklers. Farm sizes ranged from 7 to 970 ha. The average paddock size per farm ranged from 1.2 to 15.7 ha. The Northern Irriga tion District farms have the lowest average annual rainfall and the highest temperature while the West Gippsland farms have the highest average annual rainfall and a milde temperature range (Table 1). Key land management options and soil characteristics acros the farms are also summarised in Table 1. The average paddock size per farm ranged from 1.2 to 15.7 ha. The Northern Irrigation District farms have the lowest average annual rainfall and the highest temperatures, while the West Gippsland farms have the highest average annual rainfall and a milder temperature range (Table 1). Key land management options and soil characteristics across the farms are also summarised in Table 1.
Two paddocks on each of the seven calibration farms were measured regularly (weekly to monthly intervals). On each visit, three destructive samples were removed from at least five predetermined 10 × 10 m areas across the two paddocks. The 10 × 10 m areas were based on the S2 pixel grid to enable direct comparison of field data and image observations, as shown in Figure 2. The sampling procedure was replicated on the four independent validation farms. Two paddocks were sampled on each of these farms; however, these farms were sampled only once.  Three pasture observations and destructive cuts were taken in each 10 × 10 m sample pixel (sampling unit), and the mean pasture dry matter yield was calculated. The number and the orientation of the pixels across the paddocks depended on paddock size. The minimum sample per paddock was five pixels.

Ground Data Collection
Field data collection across all eleven farms was consistently applied using the following procedures. All field data collection and in-field navigation were performed using Emlid RTK GPS units with a spatial precision of +/−2 cm. The destructive sample areas Figure 2. Three pasture observations and destructive cuts were taken in each 10 × 10 m sample pixel (sampling unit), and the mean pasture dry matter yield was calculated. The number and the orientation of the pixels across the paddocks depended on paddock size. The minimum sample per paddock was five pixels.

Ground Data Collection
Field data collection across all eleven farms was consistently applied using the following procedures. All field data collection and in-field navigation were performed using Emlid RTK GPS units with a spatial precision of +/−2 cm. The destructive sample areas varied between 70 × 35 cm and 64 × 31 cm quadrats. Within these quadrats, observations of plant leaf stage, botanical composition, and three measurements of height (ruler, rising plate meter, and sonar) were taken before sampling and recorded in the field using the ESRI Collector field application [54]. The pasture was cut to ground level, bagged, weighed fresh, washed to remove soil and debris, and then oven-dried at 100 • C for 24 h or until a constant weight was recorded. In addition, a vehicle-mounted wide-angle ultrasonic sensor and a UAV-mounted multispectral camera were simultaneously used during data collection, with results published separately [55][56][57]. In the previous works, while different sensors and numerical techniques were implemented, the focus remained on the prediction of pasture yield and nutritive characteristics with similar sample pre-processing approaches. The work of Karunaratne et al. [55] focused on UAV-based empirical numerical modelling approaches by considering Structure from Motion only (Sf M), vegetation indices (VI) only, and finally, Sf M + VI combination. Meanwhile, Lawson et al. [56] presented global ultrasonic (vehicle-mounted) and RPM models using linear mixed models to predict pasture biomass, and Thomson et al. [57] compared different modelling approaches using hyperspectral datasets to predict pasture biomass and nutritive characteristics.
Data collected in the field was later downloaded from ESRI Collector via ArcGIS Online [58] (AGOL) and saved to an ESRI File Geodatabase. Polygon data were transformed from the WGS84 projection to the UTM zone appropriate for each farm to match the S2 data.

Satellite Data Collection
This study was aligned with a broader study where the aim was to test a variety of sensors. One study design and collection protocol were utilised for multiple sensors to reduce the operational cost across the programme. Since ground data collection timing depended on appropriate UAV flying conditions (wind speed and rain) and rarely coincided with satellite overpass dates, a process was developed to select the most appropriate S2 image for each collection date. This process considered images before and after field data collection, the time of year, location, grazing, and the number of days between the field data collection and image acquisition.
Intensively managed dairy pastures are grazed frequently. However, grazing times recorded by the farmer were not always available or accurate. Therefore, a method was developed to identify the timing of field measurement relative to grazing, which enabled the appropriate selection of a satellite image. The image selection technique calculated the measured mean above-ground biomass (kg DM/ha) for each paddock on each collection day. These individual collection averages were combined to create summary statistics (minimum, maximum, mean, and standard deviation) for each paddock over the complete time series. If the mean biomass on the collection day was less than the minimum plus the standard deviation of the total paddock collection biomass, then grazing had recently occurred. As grazing may have occurred the previous day, it was important to use an image acquired after the ground collection date. On the contrary, if the mean biomass on the day of collection was greater than the maximum minus the standard deviation of the total paddock collection, then it was likely that grazing would occur after the day of ground collection. Therefore, an image acquired before the ground collection date was chosen ( Figure 3). This process was only used for the calibration and internal validation dataset. Independent validation data were unique; therefore, the closest image to the collection date was chosen. Trampling the area within and around the pixel may be possible in 34% of images acquired after the ground collection date; however, this was not considered a parameter in the image selection approach. paddock collection, then it was likely that grazing would occur after the day of ground collection. Therefore, an image acquired before the ground collection date was chosen ( Figure 3). This process was only used for the calibration and internal validation dataset. Independent validation data were unique; therefore, the closest image to the collection date was chosen. Trampling the area within and around the pixel may be possible in 34% of images acquired after the ground collection date; however, this was not considered a parameter in the image selection approach.  Pasture accumulation rates differ broadly between regions and months and depend on farm management. Based on the reported accumulation rates from Figures 8, 11, 12, and 13 of Doyle et al. [59], a maximum number of days of lag per month was established to limit the potential growth differential between measurement and overpass of less than 200 kg DM/ha. The current study extended no image lag date beyond five days. Currently, no average monthly pasture growth rates are reported for irrigated dairy pastures in southeast South Australia. The southeast South Australia site was a summer irrigated system, as were the Northern Irrigation Region Sites and the Macalister Irrigation Region sites. Therefore, the shortest lag rule from those two sites was applied as a precautionary approach.
Data were collected from May 2019 to July 2021. A total of 200 paddock data collection dates and 100 S2 image datasets, with less than 15% cloud, were available across the seven calibration farms between May 2019 and March 2020. Of the available calibration data, 75 paddock collection datasets matched 40 S2 images; however, on closer inspection, the cloud affected eight images. Therefore, the final calibration data collection included 45 paddock datasets and 30 S2 images. Using the same approach for all data available from April 2020 to July 2021, a final validation collection included 16 paddock datasets and eight S2 images. Figure 4 overlays the image acquisition and ground collection dates for the calibration and validation dataset, separated by farm, for a comprehensive understanding. The reduction in data availability from March 2020 to April 2021 was due to COVID-19 lockdowns restricting farm access. the cloud affected eight images. Therefore, the final calibration data collection included 45 paddock datasets and 30 S2 images. Using the same approach for all data available from April 2020 to July 2021, a final validation collection included 16 paddock datasets and eight S2 images. Figure 4 overlays the image acquisition and ground collection dates for the calibration and validation dataset, separated by farm, for a comprehensive understanding. The reduction in data availability from March 2020 to April 2021 was due to COVID-19 lockdowns restricting farm access.

Data Processing-Satellite Spectral Index Calculations and Stacking with Bands
Level 2A, Sentinel-2A, and 2B images with less than 15% cloud were downloaded from the Sentinel Australasia Regional Access (SARA) website (https://copernicus.nci.org.au/sara. client/#/home, accessed on 15 June 2020). Image data processing occurred in ENVI ® +IDL version 5.6 [60]. Each dataset was spatially subsampled to an area of interest (AOI) that included the whole farm (not just the paddocks sampled), each band resampled to 10 m resolution and stacked in a 12-band TIF file. Information on the original S2 bands has been included in Table 2 (rows 55 to 66). The "GainOffset" function was used with offset value for all bands set to 0.0001, and a high (1.0) and low (0.0) clip was applied. Fifty-four spectral indices (SI) were then calculated for each image using an ENVI ® +IDL model and stacked into a 54-band TIF file. All SI are listed in Table 2 (rows one to 54). The initial 54 index selection was made from an extensive list curated from a literature review to identify those indices that could be calculated using S2 bands. Each SI was reviewed considering the bandwidths and centre wavelength of the S2 sensors rather than the traditional band designation or the original S2 band spatial resolution. Consideration was also given to the sensor initially used to develop the index. Each index calculation was tested in ENVI to ensure the notation was calculated as defined in the original reference. Simple ratios or difference indices were calculated using all available band combinations to test the sensitivity of the full array of bands available on the S2 sensors. Table 2. Information on the spectral indices reviewed in rows one to 54, E = Exploratory band combination, V = variation on a published spectral index. Information on original S2 bands used for the prediction model [61], in rows 55 to 66, G = Green, RE = Red Edge, NIR = Near Infrared, R = Red, BL = Blue. The final image stack index list of 21 bands is provided in column two.

Row No.
Image Index Name/Description SI Formulae or Original Band Information Source (If Applicable) Modified Chlorophyll Abs Ratio IMPROVED

Data Merging
Due to the precise locational attribution of the in-field data collection, it was possible to combine field measurements with the satellite observations into one dataset. Each ground data collection dataset with a matching S2 image was converted to a point layer using the 'Feature to Point' tool in ESRI ArcMap 10.8 [91], and a unique S2 pixel identifier was added to generate summary statistics and other analyses at pixel scale as required. Pixel-level information from the 12-band and 54-SI stack were spatially joined to the ground data point layer using the 'Extract Multi values to points' tool in ESRI ArcMap 10.8 [91]. All individual datasets were stacked using the R programming language for analysis [92]. The R computation was conducted on a Windows ® 10 operating system. The computer had a processor of 11th Gen Intel ® CoreTM i9-11950H@2.60GHz 2.61 GHz processor with 64 GB RAM. The ML was run through the NVIDIA ® RTX A3000 GPU to accelerate processing and improve numerical efficiency.

Pre-Processing and Development of Machine Learning Framework
Additional categorical variables were added to the dataset prior to modelling. As the season is an integral component of biomass modelling, the data were categorised into five seasonal periods consistent with the Australian forage value index [65,93]. The seasons were winter (June, July), early spring (August, September), late spring (October, November), summer (December, January, February), and autumn (March, April, May). In addition, information on regions and types of management (dryland or irrigated) was included (Table 1).
An ML approach was used to accommodate the high-dimensional datasets characterised by multi-collinearity among the dependent variables. The RF model was chosen among different ML approaches as it is an efficient technique for analysing linear and non-linear relationships [44]. Furthermore, RF is an ensemble, non-parametric, supervised approach that is widely applied in agricultural applications [94][95][96] and, therefore, considered appropriate to predict pasture biomass. RF models require additional functions to increase the efficiency and reliability of predictions, such as performing variable importance (VI) to optimise the feature space [44,[97][98][99], presenting a correlation matrix to observe the relationships among the high dimensional datasets [100], and the detection of outliers [101], to name a few. In addition, several statistical interpretations and analyses, along with the data preparation, were performed to improve model performance, and these have been described in the following sections.

Shapiro-Wilk Test
In ML algorithms, outliers can deceive the sample training process, leading to increased training times and less accurate models [102]. Therefore, it was essential to conduct a detection check to identify outliers and normal distribution. The Shapiro-Wilk (SW) test was considered in this study with both calibration and validation datasets to test the null hypothesis, normality, and the detection of outliers due to its proven capability in digital agriculture, mainly when dealing with multivariate datasets [103]. The null hypothesis means that no direct relationship exists between the two variables being studied. A p-value, also known as the probability value, is an important parameter that is mostly observed to accept or reject a null hypothesis. Generally, a lower p-value refers to statistical significance and if the null hypothesis is rejected. In the present study, p < 0.05 was assigned to determine statistical significance. The calculation of the SW normality test was denoted by W (0 ≤ W ≤ 1), which shows how efficiently the ordered and standardised sample quantiles fit the standard normal quantiles, with 1 demonstrating a perfect match [102]. In the present study, W = 0.9596 and W = 0.9571 were obtained for internal and independent model validations, respectively, along with p < 0.05. Therefore, it improves confidence in the development of the model.

Conditional Latin Hypercube Sampling
Conditional Latin Hypercube sampling (cLHS) can be used for the spatial splitting of datasets for model calibration and internal validation [104][105][106] as ancillary variables occupy a hypercube in the feature, attribute, or variable space. The feature space can be sampled instead of sampling the whole geographical space. However, planning a sampling methodology that can select the sampling locations covering the hypercube of the feature space is a complicated task. Therefore, a widely used alternative is to consider cLHS, which is non-parametric, and additional conditions can be considered for model optimisations based on input features [105].
The obtained data was split into calibration and validation datasets (80% [n = 171] and 20% [n = 43]) using the cLHS algorithm. The library package "clhs" was used to perform the split [105]. The evolution of the objective function for the number of iterations was monitored to observe stability. The algorithm identifies the points that can represent a Latin hypercube through a random data sampling procedure. The dataset for cLHS was subjected to 2500 iterations to reach a steady state of the objective function. Subsequently, cLHS was used as an objective function for data splitting before calibration and internal validation of the RF ML algorithm.

Variable Importance Section through Boruta Algorithm
The RF approach is efficient but could be time-consuming for application across large spatial extents. Selecting variable importance (VI) is a standard practice in ML modelling to reduce the number of input variables. At the same time, this process allows the model to emphasise the variables which are more influential in predicting the output. Since input variables cannot be reduced or eliminated arbitrarily, standard algorithms are implemented into ML models to perform the feature selection task. While there is no universal yardstick of the most accurate algorithm to perform this task, the Boruta algorithm is popular, particularly for remote sensing ML modelling [44,[107][108][109]. The Boruta algorithm is an ensemble technique and was considered in the ML approach to increase RF efficiency by reducing the number of variables. The library package "Boruta" was employed in the prediction model to perform feature selection by means of variable importance (VI) [109]. Therefore, Boruta ensures that the RF model is built based on only relevant and essential features with minimised data dimensionality.
The original variable list comprised 12 original S2 bands and 54 SI totalling 66 bands and indices ( Table 2). In addition, season and management were included in the model as categorical variables. Therefore, a total of 68 input variables were considered initially before the final model development. Multiple Boruta algorithms were performed to identify the most influential bands and SI to reduce the number of variables. The best 21 bands consisting of 12 original bands and 9 SI were unbiasedly selected, as shown in Table 2 and as described in Section 2.7. Season and management were also found to be important through the feature selection approach. Therefore, the model was developed with season and management totalling 23 input variables as selected by Boruta. Boruta considered all the candidate features (original bands and SI, season, and management) to predict pasture biomass and randomly designed shadow versions of each feature. The maximum number of iterations was varied to ensure that all essential attributes were identified and the convergence state was monitored. Boruta runs were based on iterations written as a function called "maxRuns" in R, indicating the maximum number of iterations. In the current RF model, maxRuns = 1500 was found to be adequate for identifying all important variables since the VI selection decision was reached within 1500 iterations. The RF model was optimised once all requirements were fulfilled.

Random Forest Modelling
In this research, the RF model was performed in the R statistical programming environment, and the "ranger" package was used to create the model. The "ranger" package is a fast implementation of RF modelling and is specifically well-suited for analysing high-dimensional data [110,111]. Furthermore, this RF model was simple to implement, and in the current study, two key model parameters were needed to define: Ntree and mtry [110], along with additional associates to ensure smooth execution. The following briefly discusses some of the important parametric values considered in the RF model: RF model parameters: The primary objective of RF model parameter tuning is to find the optimal hyperparameters for the calibration dataset. While the default RF setting has been reported to be less tunable, specific tuning search strategies (such as grid search) have assisted in evaluating the discrete parameters for the out-of-bag performance [110]. One of the primary tools is the "train" function embedded with the "caret" package [112], which was used to evaluate the RF model parameter tuning and estimate the performance from the training set. Furthermore, "tuneGrid" and "trainControl" arguments were integrated into the coding environment to generate the parametric values of the candidates and modify the resampling (if required). Additionally, the "repeatedcv" argument was used to integrate the repeated K-fold cross-validation, where the argument repeats controlled the number of repetitions. A specific number maintains K in the argument; by default, it is 10.
• Ntree: The parameter Ntree refers to the number of decision trees generated. As per the literature, the standard requirement to analyse remote sensing data, i.e., a default value of Ntree = 500 was used [44]. • mtry: The parameter mtry denotes the number of variables to be selected and tested for the best split while growing trees. Lower mtry values have been attributed to stability enhancement, as it reduces the number of correlated trees [110]. Several tests were included before selecting mtry values, as advised by Probst et al. [110]. In the present model, mtry = 6 was found to be the best with reduced computational time.
In RF, the default splitting rule considers a selection of variables with many possible combinations. However, it has been reported that the default rule favours the selection of the variables with many splits and can overlook variables with fewer splits [110]. This type of variable bias decreases RF efficiency, so a modified approach has been considered in this study. In this study, the splitting rule was randomised per the technique described in Geurts and Wehenkel [113]. In addition, the "extratrees" feature was considered. The function "extratrees" adds additional components of randomness to the trees in RF modelling. However, it has been observed that such randomness substantially impacts VI ranking and may lead to RF being biased [110]. The "permutation" function of VI was integrated into the R code to ensure an impartial VI [114].
A complete flowchart outlining the data manipulation and analysis framework using the ML approach is shown in Figure 5.
trolled the number of repetitions. A specific number maintains K in the argument; by default, it is 10.
• Ntree: The parameter Ntree refers to the number of decision trees generated. As per the literature, the standard requirement to analyse remote sensing data, i.e., a default value of Ntree = 500 was used [44]. • mtry: The parameter mtry denotes the number of variables to be selected and tested for the best split while growing trees. Lower mtry values have been attributed to stability enhancement, as it reduces the number of correlated trees [110]. Several tests were included before selecting mtry values, as advised by Probst et al. [110]. In the present model, mtry = 6 was found to be the best with reduced computational time.
In RF, the default splitting rule considers a selection of variables with many possible combinations. However, it has been reported that the default rule favours the selection of the variables with many splits and can overlook variables with fewer splits [110]. This type of variable bias decreases RF efficiency, so a modified approach has been considered in this study. In this study, the splitting rule was randomised per the technique described in Geurts and Wehenkel [113]. In addition, the "extratrees" feature was considered. The function "extratrees" adds additional components of randomness to the trees in RF modelling. However, it has been observed that such randomness substantially impacts VI ranking and may lead to RF being biased [110]. The "permutation" function of VI was integrated into the R code to ensure an impartial VI [114].
A complete flowchart outlining the data manipulation and analysis framework using the ML approach is shown in Figure 5.

Spectral Index Reduction
After several primary iterations for biomass estimation and prediction model development, 12 bands and 9 SI were selected using the outcomes of the Boruta algorithm. At least 12 iterations were made, and the results of the variable importance plot were collated. Row numbers 3, 14, 16, 22, 32, 34, 43, and 54 (column 1, Table 2) were always rejected. Row numbers 6,8,12,13,17,19,23,36, and 48 (column 1, Table 2) were always important and used in the biomass estimation and prediction model. The remaining 37 spectral indices had variable results across model iterations. All were rejected at least once, and some up to eight times.

Model Quality Assessment
As mentioned earlier, the model prediction quality was tested through internal and fully independent validation. Different metrics were used to evaluate the model quality: (a) RMSE, which provides information on the model accuracy by calculating the differences between the predicted and observed values, (b) NRMSE, which is the calculated RMSE as a percentage of the measured mean of the data, (c) R 2 , and finally (d) Lin's concordance correlation coefficient (LCCC).
The following equations describe the expression to determine the validation indices: where Y = measured pasture DM yield, Y i = predicted pasture DM yield, and N is the number of samples. NRMSE = (RMSE/Mean measured pasture DM yield) × 100 (2) LCCC depicts the variation of the predicted values from a unity line (1:1) [115]. LCCC is a measurement of the degree to which the predicted values adhere to the concordance line of slope 1.0 through the origin and an outcome of the product between the Pearson correlation coefficient and a bias factor reflecting the precision and accuracy. The bias factor information is achieved from mean and slope bias [115]. LCCC is defined as the following: where m x and m y are the measured and predicted pasture biomass means, respectively; s x 2 and s y 2 represent the variances of the measured and predicted pasture biomass. Meanwhile, r is the Pearson correlation coefficient between measured and predicted pasture biomass.

Pearson Correlation Matrix and Best-Performing ML Model
Examining the relationships among the variables as well as the feature selection are two of the significant steps in the ML approach. Importantly, weaker relationships between the predictor and target variables affect the efficiency, prolong the training times of the samples, increase the computational complexities, and reduce the potential impact of significant parameters. Pearson's correlation matrix is a popular method for demonstrating correlation among variables in ML. Each cell in the correlation matrix represents the correlation between two variables. The values are located between −1 and +1, with values closer to −1 and +1 indicating strong negative and positive correlations. A value closer to 0 demonstrates no connection between those variables. This study calculated Pearson's correlation matrix to examine the correlation coefficient among the selected bands and spectral indices, as shown in Figure 6. In general, it was observed that all the bands considered to predict mean values of the dry matter showed good correlations, with Image Index Band_20 (B11 SWIR), Band_11 (B2 Blue), Band_19 (B9 Water vapour), and Band_21 (B12 SWIR) showing the strongest correlations. None of the bands exhibited a null value (0) to predict dry matter. Figure 7a shows the internal model validation where R 2 = 0.90 and LCCC = 0.72 were obtained. Moreover, independent validation (n = 84) was conducted to prove the model's ability to be a global model. Figure 7b shows the correlation between the measured and predicted values, where R 2 = 0.88 and LCCC = 0.68 were obtained. The R 2 and LCCC values suggested that the model could adequately predict a fully independent validation dataset. correlation matrix to examine the correlation coefficient among the selected bands and spectral indices, as shown in Figure 6. In general, it was observed that all the bands considered to predict mean values of the dry matter showed good correlations, with Image Index Band_20 (B11 SWIR), Band_11 (B2 Blue), Band_19 (B9 Water vapour), and Band_21 (B12 SWIR) showing the strongest correlations. None of the bands exhibited a null value (0) to predict dry matter.  Table 2, column 2). DM stands for dry matter. A dark thin diagonal straight line represents a perfect 1:1 relationship, and as the relationship reduces, the line becomes wider and lighter in colour. Figure 7a shows the internal model validation where R 2 = 0.90 and LCCC = 0.72 were obtained. Moreover, independent validation (n = 84) was conducted to prove the model's ability to be a global model. Figure 7b shows the correlation between the measured and predicted values, where R 2 = 0.88 and LCCC = 0.68 were obtained. The R 2 and LCCC values suggested that the model could adequately predict a fully independent validation dataset.

Combined Validation
Pasture biomass measured in all seasons and types of management (dryland/irrigated) was considered collectively in both validation datasets. The maximum yield range of the internal validation dataset was close to 6000 kg DM/ha (Figure 7a), which was larger Pasture biomass measured in all seasons and types of management (dryland/irrigated) was considered collectively in both validation datasets. The maximum yield range of the internal validation dataset was close to 6000 kg DM/ha (Figure 7a), which was larger than the independent validation dataset. Although the maximum predicted value was not as close to the regression line, there were limited data at this high range. Furthermore, the prediction range was diminished in the independent validation due to limited relevant data (locations and seasons) (Figure 7b). Nevertheless, the maximum range of biomass prediction exceeded 3000 kg DM/ha in both internal and independent validation (Figure 7), which provides confidence in the robustness of the model. The values of R 2 and LCCC for both validations indicate the versatility of the model to be one single global model. than the independent validation dataset. Although the maximum predicted value was not as close to the regression line, there were limited data at this high range. Furthermore, the prediction range was diminished in the independent validation due to limited relevant data (locations and seasons) (Figure 7b). Nevertheless, the maximum range of biomass prediction exceeded 3000 kg DM/ha in both internal and independent validation ( Figure  7), which provides confidence in the robustness of the model. The values of R 2 and LCCC for both validations indicate the versatility of the model to be one single global model.   Quantitative comparisons are shown in Table 3. The minimum and maximum measured biomass are shown to highlight the extensive range of biomass values considered in the model. The w and p values from the SW test also show that the data considered for the model calibrations were statistically significant. Model accuracy indicators, RMSE and NRMSE, provide additional information on the robustness of the model. An RMSE < 500 kg DM/ha was achieved overall (RMSE = 439.49 kg DM/ha and 457.05 kg DM/ha for internal and independent validations, respectively). The entire independent validation data were entered directly into the already calibrated model to assess accuracy, as shown in Table 3.

Validations Based on Season and Management
The VI plot ( Figure 8) indicated that season was the most important variable. However, it should be mentioned that while the ranks of variables differed, all the variables in Figure 8 were important in the biomass prediction model and should not be eliminated. The type of farm management (irrigated or dryland) was another important variable that passed the Boruta feature selection and was important for the biomass prediction model. Therefore, separate plots were generated for all seasons and farm management (dryland or irrigated) to assess the model accuracy separately. Quantitative comparisons are shown in Table 3. The minimum and maximum measured biomass are shown to highlight the extensive range of biomass values considered in the model. The w and p values from the SW test also show that the data considered for the model calibrations were statistically significant. Model accuracy indicators, RMSE and NRMSE, provide additional information on the robustness of the model. An RMSE < 500 kg DM/ha was achieved overall (RMSE = 439.49 kg DM/ha and 457.05 kg DM/ha for internal and independent validations, respectively). The entire independent validation data were entered directly into the already calibrated model to assess accuracy, as shown in Table 3.

Validations Based on Season and Management
The VI plot ( Figure 8) indicated that season was the most important variable. However, it should be mentioned that while the ranks of variables differed, all the variables in Figure 8 were important in the biomass prediction model and should not be eliminated. The type of farm management (irrigated or dryland) was another important variable that passed the Boruta feature selection and was important for the biomass prediction model. Therefore, separate plots were generated for all seasons and farm management (dryland or irrigated) to assess the model accuracy separately.  Table 2 for the description of Bands. Season and management are categorical variables and are described in Section 2.6.3. Figure 8. The variable importance generated from the best-performing RF model through the Boruta algorithm. Refer to Table 2 for the description of Bands. Season and management are categorical variables and are described in Section 2.6.3.

SWIR Band Validation
Several tests across different farms were performed to determine the importance of SWIR bands in pasture biomass estimation. SWIR bands were partially or completely removed from the variable set through many model iterations. Index Number Band_20 (B11 SWIR), Band_21 (B12 SWIR), and Band_3 (B12 SWIR/B8A NIR) were left out of the test to determine the importance of SWIR bands. The internal validation results are presented in Figure 11. The removal of SWIR bands reduced the accuracy of the model significantly. The accuracy indicators show that R 2 = 0.79 and LCCC = 0.57 were obtained, which are approximately 12% and 21% reductions from the original prediction model. In addition, RMSE = 635.46 kg DM/ha was achieved through this test, which was approximately 45% augmentation on the original pasture biomass prediction model presented in Section 3.2 and Table 4.

SWIR Band Validation
Several tests across different farms were performed to determine the importance of SWIR bands in pasture biomass estimation. SWIR bands were partially or completely removed from the variable set through many model iterations. Index Number Band_20 (B11 SWIR), Band_21 (B12 SWIR), and Band_3 (B12 SWIR/B8A NIR) were left out of the test to determine the importance of SWIR bands. The internal validation results are presented in Figure 11. The removal of SWIR bands reduced the accuracy of the model significantly. The accuracy indicators show that R 2 = 0.79 and LCCC = 0.57 were obtained, which are approximately 12% and 21% reductions from the original prediction model. In addition, RMSE = 635.46 kg DM/ha was achieved through this test, which was approximately 45% augmentation on the original pasture biomass prediction model presented in Section 3.2 and Table 4.
A comparison table has also been presented in Table 4 to demonstrate the importance of SWIR bands by comparing two models, one with SWIR bands and one without. The model with SWIR bands predicted maximum pasture biomass up to 4348.25 kg DM/ha from a specific day in the summer season of the PS04 farm. The alternative model, excluding SWIR bands (Band_20, Band_21) and Band_3, was used to predict the pasture biomass of the same day and location, and the maximum biomass obtained was 3379.62 kg DM/ha, which was approximately 22% less than the original prediction model.  A comparison table has also been presented in Table 4 to demonstrate the importance of SWIR bands by comparing two models, one with SWIR bands and one without. The model with SWIR bands predicted maximum pasture biomass up to 4348.25 kg DM/ha from a specific day in the summer season of the PS04 farm. The alternative model, excluding SWIR bands (Band_20, Band_21) and Band_3, was used to predict the pasture biomass of the same day and location, and the maximum biomass obtained was 3379.62 kg DM/ha, which was approximately 22% less than the original prediction model.

Model Automation
The model introduced in this study was used in an automated data processing flow. Farms involved in the study were mapped, and each paddock management was noted. It was common for dryland farms to have a small proportion of irrigated paddocks and vice versa. The SARA database (https://copernicus.nci.org.au/sara.client/#/home, accessed on 6 October 2020) was reviewed daily using the farm boundary, and when an image with less than 15% cloud was identified, the image was downloaded. This review and download process was automated using Python and APIs established by Geoscience Australia. Once downloaded to an internal computer, the images were pre-processed using the ESA SNAP© graph processing tool. Some of the major S2 image processing operations involved resampling, clipping, S2 band arithmetic computation, band merging, and finally, writing in TIF format.
The image was clipped to the farm extent, all bands were resampled to 10 m, the nine important SIs were calculated, and a stack of 12 bands and 9 SI was created. The model was augmented with functions that read the season from the image file name and the management of each paddock from the farm mapping geodatabase. Finally, biomass predictions were made per pixel to accommodate the different management types in each paddock. Prediction per pixel was achieved using the "splancs" library in the R programming language, which performed the spatial point-pattern analysis [116]. Furthermore, the "inpip" function inside the model selected points inside a polygon based on an area of interest and geospatial information. The output prediction map integrated dryland and Figure 11. Prediction model results without including SWIR bands and associated SI.

Model Automation
The model introduced in this study was used in an automated data processing flow. Farms involved in the study were mapped, and each paddock management was noted. It was common for dryland farms to have a small proportion of irrigated paddocks and vice versa. The SARA database (https://copernicus.nci.org.au/sara.client/#/home, accessed on 6 October 2020) was reviewed daily using the farm boundary, and when an image with less than 15% cloud was identified, the image was downloaded. This review and download process was automated using Python and APIs established by Geoscience Australia. Once downloaded to an internal computer, the images were pre-processed using the ESA SNAP© graph processing tool. Some of the major S2 image processing operations involved resampling, clipping, S2 band arithmetic computation, band merging, and finally, writing in TIF format.
The image was clipped to the farm extent, all bands were resampled to 10 m, the nine important SIs were calculated, and a stack of 12 bands and 9 SI was created. The model was augmented with functions that read the season from the image file name and the management of each paddock from the farm mapping geodatabase. Finally, biomass predictions were made per pixel to accommodate the different management types in each paddock. Prediction per pixel was achieved using the "splancs" library in the R programming language, which performed the spatial point-pattern analysis [116]. Furthermore, the "inpip" function inside the model selected points inside a polygon based on an area of interest and geospatial information. The output prediction map integrated dryland and irrigated biomass predictions in one farm prediction. A flowchart is presented in Figure 12, outlining the steps.
The automation process did not require user-defined commands after a farm was mapped and included in a farm geodatabase (ESRI). An example of the results across five different seasons for the PS04 farm is shown in Figure 13. It could be observed that the maximum biomass prediction exceeded 4000 kg DM/ha for early spring (Figure 13b), late spring (Figure 13c), and summer ( Figure 13d). However, the maximum range was found to be 3633.83 kg DM/ha in the winter season ( Figure 13a) and 3698.55 kg DM/ha in autumn (Figure 13e). This example was from an irrigated farm. Clear differences in biomass prediction can be seen throughout the image. Many are associated with farm infrastructure, such as water troughs, channel banks, and fences delineating paddock boundaries. Importantly, within-paddock variation is also visible due to uneven pasture growth. In this farm, uneven pasture growth is due to water availability rather than differences in soil type. This farm is flood-irrigated, and therefore, most paddocks have been levelled to allow for optimal water flow. On occasion, water does not reach the end or edges of the paddock due to low water levels in the channel, causing insufficient flow down the paddock. Due to the price and availability of water, not all paddocks were irrigated during this time, highlighting the importance of accommodating changes in management and the full range of biomass amount. irrigated biomass predictions in one farm prediction. A flowchart is presented in Figure  12, outlining the steps. The automation process did not require user-defined commands after a farm was mapped and included in a farm geodatabase (ESRI). An example of the results across five different seasons for the PS04 farm is shown in Figure 13. It could be observed that the maximum biomass prediction exceeded 4000 kg DM/ha for early spring (Figure 13b), late spring (Figure 13c), and summer ( Figure 13d). However, the maximum range was found to be 3633.83 kg DM/ha in the winter season ( Figure 13a) and 3698.55 kg DM/ha in autumn ( Figure 13e). This example was from an irrigated farm. Clear differences in biomass prediction can be seen throughout the image. Many are associated with farm infrastructure, such as water troughs, channel banks, and fences delineating paddock boundaries. Importantly, within-paddock variation is also visible due to uneven pasture growth. In this farm, uneven pasture growth is due to water availability rather than differences in soil type. This farm is flood-irrigated, and therefore, most paddocks have been levelled to allow for optimal water flow. On occasion, water does not reach the end or edges of the paddock due to low water levels in the channel, causing insufficient flow down the paddock. Due to the price and availability of water, not all paddocks were irrigated during this time, highlighting the importance of accommodating changes in management and the full range of biomass amount.

Discussion
This study has presented results that are new to the research area. The study has combined a complete dataset collected from active commercial dairy farms with a modelling approach that emphasises the crucial role of remote sensing in solving one of the significant challenges faced by the pasture-based dairy industry, i.e., automated biomass

Discussion
This study has presented results that are new to the research area. The study has combined a complete dataset collected from active commercial dairy farms with a modelling approach that emphasises the crucial role of remote sensing in solving one of the significant challenges faced by the pasture-based dairy industry, i.e., automated biomass prediction with high accuracy for dense, moist, fast-growing canopies of pasture forages. Excellent prediction qualities have been achieved for all regions, seasons, and management types, producing a single global model for automated pasture biomass prediction. Whilst other studies have investigated minor components common to this study, none have demonstrated real-world outcomes or on a scale that can definitively identify the critical spectral bands required for an accurate prediction. The following sections discuss these novel outcomes in detail.

Overview of the Prediction Model Accuracy
This study has successfully achieved the stated aim of predicting perennial ryegrass biomass with an accuracy of +/−500 kg DM/ha using internal (RMSE 439.49 kg DM/ha) and independent validation (RMSE 457.05 kg DM/ha). In addition, separate prediction accuracies were presented for seasonal and farm management differences to validate the model's accuracy. Internal validation across all seasons suggests that the model obtained R 2 ≥ 0.81 and LCCC ≥ 0.59. Independent validation was conducted across three seasons, with the model showing an accuracy of R 2 ≥ 0.83 and LCCC ≥ 0.66. Similar validations for dryland and irrigated farm management systems provided further evidence of the robustness of the present model where R 2 ≥ 0.86 and LCCC ≥ 0.67 (internal validation) and R 2 ≥ 0.90 and LCCC ≥ 0.72 (independent validation) were obtained. It should be noted that the validation data was predominantly acquired from entirely different farms and at least a year after the calibration data collection. Additionally, all data were sourced from commercial farms in paddocks under grazing, based on the specific farmer practices. Therefore, there was almost always a day offset between data collection and image acquisition. Nevertheless, R 2 and LCCC showed excellent agreements between the predicted and measured pasture biomass for internal and independent validation despite the variability of the conditions and environments under which data were collected.

Significance of S2 SWIR Bands in Improving Prediction Accuracy
The SWIR validation results were included in Figure 11 and Table 4 to demonstrate the significance of SWIR bands. They support the hypothesis that using SWIR bands would improve the ability to predict pasture yields above 3000 kg DM/ha and, therefore, enable measurements of high-yielding pastures at any stage in the growth cycle in irrigated and dryland production. Dairy pastures before harvest are dense, high moisture swards that can often exceed LAI > 3 and quickly reach the saturation point for common vegetation indices [41]. Saturation is evident at levels >0.7 NDVI [117] and biomass saturation at 2500 kg DM ha [11]. While Edirisinghe et al. [30] and Sinde-Gonzales et al. [12] reported calibration data that exceeded 3000 kg DM/ha, prediction ranges rarely exceeded these levels using NDVI-like indices. SWIR bands are known to be influenced by plant water concentration and therefore are likely to extend the saturation sensitivity point [118,119]. The importance of SWIR bands in improving biomass prediction was also supported by Dang et al. [120], Numata et al. [36], and Pandit et al. [121]. Furthermore, SWIR bands were reported to have strong correlations with pasture biomass, regardless of the type and density of the data and environmental factors [122][123][124]. This study found that adding SWIR bands increased the predictable range of biomass above 4000 kg/DM ha.
Of the 23 bands identified in the VI plot process, S2 B11 SWIR and S2 B12 SWIR were ranked second and fifth, respectively, further supporting the hypothesis that the utilisation of SWIR bands would be significant for model prediction. These bands have centre wavelengths of approximately 1600 and 2200 nm, respectively. The 1600 nm range relates to biochemical components, starch, and sugar, and the 2200 nm range to structural components, lignin, and cellulose [125,126]. Among other bands, in the top 10 most important variables, S2 B2 Blue and S2 B9 Water Vapour were present. S2 B2 has a centre wavelength of 490 nm, an area of strong chlorophyll b absorption [127]. S2 B9 has a centre wavelength of 940 nm and measures water vapour's absorption over land [128]. Both are important areas of sensitivity in dense, moist canopies. Importantly, this study found that the five most important variables did not include the commonly used red or near-infrared bands found in many biomass studies. Variables ranked six to ten had these common bands but as part of an SI rather than the original band; moreover, the S2 B8 NIR is the lowest-ranked variable in the VI plot. The band ranking suggested that bands associated with canopy and plant moisture were the most important for predicting dense, moist, high-biomass vegetation, such as PRG-based dairy pastures.

Consideration of ML for Data Analysis and Model Development
This study investigated an extensive range of spectral indices and satellite bands using the ML RF model technique to predict pasture biomass. The ML RF model accommodates high-dimensional data [44]. It was essential to use an approach that could accommodate many variables (68 in the original model) due to the complexity of the farming environment under investigation. Conventional statistical models generally fail to meet accuracy requirements in biomass estimation, particularly at the farm scale, to support on-farm decision-making [129,130]. For example, Ali et al. [129] achieved R 2 between 0.57 and 0.86 by using ML to predict the grassland yield estimation through satellite data. In contrast, the conventional statistical model demonstrated the best R 2 of 0.31. The advantage of an ML model by Ali et al. [129] was shown further because the lowest NRMSE was reported to be 25.05 using a conventional statistical model and 11.07 using an ML model. However, it should be noted that Ali et al. [129] only considered a single farm with an area of 100 ha, and yet they could only obtain R 2 values at most 0.86. Yang et al. [48] also showed a similar comparison and obtained R 2 and RMSE of 0.85 and 355 kg DM/ha for the best-performing ML model compared to R 2 = 0.58 and RMSE = 536 kg DM/ha through the traditional statistical model. Due to the large number of variables used in the present study, a direct comparison between a conventional statistical approach and an ML approach was not possible. This limitation proves that the ML approach has enhanced the prediction accuracy of the study by effectively accommodating a large data set with 23 input variables.

Impact of Soils, Climate, and Farm Activities on Satellite Images and Biomass Estimation
This study has benefited from using on-ground data covering various soils, management practices, regions, seasons, and years. The number and type of data points have provided a wide range of biomass not found in other studies. For example, the analysis presented by Chen et al. [25] did not consider large-scale farm areas, and the pasture biomass prediction barely went beyond 4000 kg DM/ha. In addition, the model quality indicators were less accurate than the present study. Chen et al. [25] reported 0.63 ≤ R 2 ≤ 0.68 and 0.22 ≤ R 2 ≤ 0.40 through their ML internal and independent validations, respectively. Wang et al. [131] obtained R 2 = 0.67 in their best-performing ML models; however, the work did not involve any grazing during the growing season. The earlier work by Numata et al. [36] considered grazing intensity, assuming a minimum effect of grazing rotation on the pasture biophysical variables of an actual farm; however, this study does not consider this as a typical scenario for a commercial farm. Xu et al. [47] presented different regression models for different regions. Each region had three equations (linear, power, and exponential) with no possible resolution toward a global or single prediction model within the area of interest. The mathematical interpretation by Grigera et al. [37] to build a correlation between NDVI and fraction absorbed by the canopy (fPAR) was derived from the literature and parametrised based on assumptions. In addition, the determination of radiation use efficiency (RUE) was limited to the local conditions with restricted time steps, and the empirical correlations require further clarification. As such, the correlation presented by Grigera et al. [37] was limited to certain locations, the quality of the data in the literature, and specific time scales. The model will not have the potential to become a global model for predicting pasture biomass. Crabbe et al. [21] primarily relied on field measurements over October and February; therefore, the effect of seasonal cycles was still limited in their prediction model. While comparing different ML techniques was interesting, the precise information on pasture biomass and other model quality indicators, such as RMSE and LCCC, was not reported. Based on the discussions above, the present study fairly outperformed those models by considering one single prediction model across all regions with an entire annual cycle, different types of farm management, and large-scale active farm data that included grazing.
Of the data collected, internal and independent validations demonstrated insensitivity to farm management (dryland and irrigated) and environmental factors (soil type and colour). This result was somewhat surprising, considering that 11 commercial farms in this study had differing soil types and colours. The VI plot assessment (Figure 8) determined that management (dryland and irrigated) was the third least important variable (19th of 21 variables). Separate validation plots were presented in Figure 10, where good accuracy was established in both dryland and irrigated management, indicating the robustness of the prediction model for large-scale working farms, regardless of the kind of management. The soil types in the farms varied from dense and sodic subsoils to strongly acidic brown chromosols, wet soils such as hydrosols, and friable and structured iron-rich soils such as red Ferrosols. Furthermore, variations in soil colour will be exacerbated with varying soil moisture. For example, a dry brown chromosol will be different in colour compared to a wet brown chromosol. Since dryland and irrigated farm management were considered for the model development, substantial variations of soil colour have been accommodated in the final model. Even so, the model was still able to obtain high accuracy, indicating that the prediction model was not biased toward any specific type of soil (soil colour or properties) or geographic location (such as southeast Australia).

Limitations of the Model
While management was of little importance to the model, the season was the most important variable identified in the VI plot of the best-performing model. Correlations for each of the five seasons used were significant (>0.8). The lowest correlation was found for the summer season (Dec, Jan, and Feb), where the most extensive data range was measured. High biomass occurred in summer irrigated pastures while dryland pastures commonly senesce due to lack of moisture and have low biomass at this time. While Figure 9 showed clear groupings associated with each season demonstrating good model outcomes, further improvement is likely if the climatic characteristics recorded at the overpass were used rather than coarse seasons. The enhancement will have the most significant impact in the summer when the difference in soil moisture concentration and evapotranspiration is most marked between irrigated and dryland management systems.
Even though this study has demonstrated the transferability of the model across regions within southeastern Australia, seasons, and management types due to an extensive dataset, the present model was limited by cloud-free imagery. While seasons like late spring (October, November) and summer (December-February) had numerous cloud-free satellite overpasses, the major challenge was during the winter season (June-July). As a result, the range of biomass estimation was limited during the winter. Furthermore, the satellite images were only considered if the field data collection occurred after S2 overpassing on the same day or within a time (≤5 days) that limited potential pasture growth to less than 200 kg DM/ha. Additionally, the present study contains field data from commercial farms; therefore, activities such as grazing events will also affect the paddock's variability. Furthermore, grazing only sometimes follows a sequential schedule, depending on pasture growth, suitable weather, and farm requirements. Consequently, available satellite data might be further impacted. The ideal approach was to coincide the field data collection date with the same day of the satellite overpass. However, this was difficult to achieve due to certain farm activities, weather conditions, and set satellite overpass schedules. Among the whole dataset considered in this study, less than ten field data collections coincided with the satellite overpassing, which might also have affected the model prediction accuracy.
For pasture to be the cheapest feeding source for dairy animals, grazing must be well-managed to ensure efficient pasture utilisation. Therefore, the limited temporal availability of biomass estimations may significantly impact the practical application of the model outcomes developed. The addition of Synthetic-aperture radar (SAR), such as Sentinel-1 (S1), which can pass through the cloud and provide regular observation, is under development to overcome this significant issue. However, SAR satellites frequently face challenges due to speckle, which is caused by the interference between the randomly distributed scatters within each pixel [132]. Hence, building a multivariate hybrid model requires more investigation into image pre-processing, calibration, and integration with the S2 imagery. A similar ML RF modelling approach integrating S1 data will be investigated based on the model's success described here.

Conclusions
This study has presented an ML RF model incorporating S2 images to predict PRG biomass in pastures across southeast Australia on 11 commercial dairy farms. This study has successfully achieved the stated aim of predicting perennial ryegrass biomass with an accuracy of +/−500 kg DM/ha using internal and independent validation. R 2 and LCCC showed excellent agreements between the predicted and measured pasture biomass for internal and independent validation despite the variability of the conditions and environments under which data were collected.
This study has also found that the availability of SWIR bands increased the predictable range of biomass above 4000 kg/DM ha. The R 2 value and LCCC showed excellent agreement between the predicted and measured pasture biomass amount (kg DM/ha) for internal and independent validations despite the complexity of the data. The results support the hypothesis and exceed the expectation that using SWIR bands would improve the ability to predict pasture yields above 3000 kg DM/ha and, therefore, enable measurements of high-yielding pastures at any stage in the growth cycle in irrigated and dryland production. The VI plot band ranking suggested that bands associated with canopy and plant moisture were the most important for predicting dense, moist, high-biomass vegetation such as PRG-based dairy pastures.
The ML RF approach improved the prediction accuracy of the study by effectively accommodating a large set of variables (12 bands, 9 SI, season, and management). This study outperformed previously reported models by considering one single prediction model across all regions with a complete annual cycle, different types of farm management, and typical commercial farm activities, including grazing. Additionally, the model was still able to obtain high model quality, indicating that the prediction model was not biased toward any specific type of soil (soil colour or properties) or geographic location.
Further improvement is possible if the climatic characteristics recorded during the overpass were used rather than coarse seasons. The enhancement will have the most significant impact in the summer when the difference in soil moisture concentration and evapotranspiration is the greatest between irrigated and dryland management systems. Among the whole dataset considered in this study, less than ten field data collections coincided with the satellite overpassing, which might also have affected the model prediction accuracy. Based on the model's success described here, a similar ML RF modelling approach will be investigated to integrate S1 data.  Data Availability Statement: Data may be made available on request at the completion of the research program when all publication is complete.