Geographically weighted spatial modelling of sediment quality in Lake Okeechobee, Florida

Abstract

Lake Okeechobee is a large shallow lake with over 44% (as of 2006) being underlain with phosphorus-enriched mud sediments. Wind-induced sediment re-suspension may play a significant role in the nutrient cycling of this lake. It is critical to develop optimal models with low uncertainty to map the distribution and changes of total phosphorus (TP) over time. Geographically Weighted Regression (GWR) is a powerful spatial regression method that examines the details of relationship between the target variable and independent variables and their changes over space. The calibrated GWR models were applied to data sets of organic mud sediment from Lake Okeechobee collected in 1988, 1998 and 2006. Three sets of ancillary data were also used in the models: total iron (Fe), which was strongly and positively correlated to TP (0.69–0.72); mud thickness, which was moderately and positively correlated to TP (0.57–0.72); and site elevation, which was weakly but negatively correlated to TP (−0.46- −0.55). The GWR models were compared with Ordinary Kriging (OK), Ordinary Least-squares Regression (OLS) and Regression-Kriging (RK) to determine their accuracy. GWR models use both spatial auto-correlation and correlation between TP and independent variables to improve the model performance. The GWR models were superior to OK and RK models previously applied to these data sets. The GWR (total Fe) and GWR (mud thickness and site elevation) models were the most accurate based on lower root mean square errors (RMSE). GWR (Fe) and GWR (Thickness and Elevation) models were used to map TP concentrations and TP mass. Using TP concentrations and mud volumes, TP mass in the lake sediments was estimated for each sample set. The TP mass increased about 38–44% from 1988 to 1998 and decreased about 30–34% from 1998 to 2006. TP mass was similar in 1988 and 2006.

Keywords:

• geographically weighted regression (GWR)
• Lake Okeechobee
• model accuracy
• ordinary kriging (OK)
• TP mass
• total phosphorus (TP)
• spatial models

Introduction

Lake Okeechobee is a large shallow lake with over 44% of the lakebed containing phosphorus-enriched mud sediments (Yan and James 2007). Because of the strong link between water quality and resuspended sediments (James and Zhang 2008), the distribution of these sediments, their compositions and changes over time, sediment re-suspension and interaction with water quality are of direct concern for lake management and restoration efforts to improve the lake’s environmental quality. In 1988, a comprehensive survey reported that the upper 10 cm of mud sediments within the lake contained an estimated 28,600 metric tons of P (Reddy, Sheng, and Jones 1995). A second survey in 1998 produced similar estimates (Fisher, Reddy, and James 2001). Both studies indicated a large variation (coefficient of variation greater than 28%) in these estimates.

In order to improve the estimates of P content and distribution of P in the sediments, generic and robust spatial modelling techniques can be used. Kriging and its variants have been widely recognized as primary spatial interpolation techniques from the 1950s. In the 1990s, with the emergence of geographical information science (GIS) and remote-sensing technologies, exhaustively mapped ancillary variables (e.g. digital elevation terrain attributes), images and other derivatives (such as land-use data) are used to estimate soil parameters directly with linear regression models (Gessler et al. 1995; Moore et al. 1993). This approach is termed environmental correlation (McKenzie and Ryan 1999) since quantitative environmental variables were used in the regression. In the last decade, a number of ‘hybrid’ interpolation techniques, which combine kriging and ancillary variables, have been developed and tested. Kriging combined with regression (McBratney et al. 2000) is considered as one of the most popular methods with few estimated model parameters (Knotters, Brusm, and Oude Voshaar 1995; Hengl 2009). Most recently, an optimized regional regression (ORR) model was proposed to estimate house value and it produced the most accurate predictions compared to four other regression methods (Lu et al. 2013).

Pebesma (2004) selected models for spatial prediction based on three questions: (1) is a physical model defined; (2) do sampled variables correlate with environmental variables and (3) do residuals show spatial autocorrelation? Four potential models can be chosen depending on how the questions are answered for the area of interest (Figure 1). If a deterministic model is not defined, and correlations exist between sampled variables and environmental factors, then a statistically significant regression model is defined. The regression residuals are then examined for spatial autocorrelation. If no spatial autocorrelation is identified with the residuals, the regression coefficients can be estimated using ordinary least squares (OLS). Otherwise, mixed models like regression-kriging (RK) and Geographically Weighted Regression (GWR) can be applied. If the data are not correlated to environmental factors but do exhibit spatial autocorrelation, the variogram of the target variable (describing the degree of spatial dependence of the variable) can be analysed and OK and its variants can be used for estimation. Otherwise, a global mean is the best estimate for the area of interest.

Figure 1. Model selection processes for spatial prediction (OLS: ordinary least square, OK: ordinary kriging; RK: regression-kriging and GWR: geographically weighted regression).

The objective of this study is to develop optimal spatial models for accurate estimates of total phosphorus (TP) concentrations and mass in Lake Okeechobee sediments over time (1988–2006). Spatial changes also are examined to identify the potential effects of extreme environmental forcing during 2004 and 2005 hurricane seasons. Several unique GWR models were calibrated and validated to describe the sediment characteristics, and then the best models were used to calculate TP and the changes in TP from 1988 to 2006.

Materials and methodology

Study area and data

Lake Okeechobee is a large shallow lake with an area of approximately 1730 km² and an average depth of 2.7 m (Figure 2). It is the largest freshwater lake in the southern United States. It is located in the southern peninsula of Florida with a subtropical climate. Average rainfall is 1.12 m (SFWMD 2014), evaporation is in excess of 1.3 m (Abtew 2001) and water temperatures range from 10°C on rare occasions in winter to 32°C in summer. The lake experiences major storm events once or twice every decade (James and Pollman 2011).

Figure 2. Bathymetry, zones and sample sites of Lake Okeechobee.

Lake Okeechobee is turbid and has experienced decades of excessive phosphorus loading, primarily from agricultural runoff. This excessive P load has accumulated primarily in mud sediments that comprise 44% of the lakebed (James and Pollman 2011). While three major regions of the lake are easily defined – the pelagic, nearshore and littoral (Figure 2) – there are distinct regions comprising different types of sediments. In addition to the mud sediments, three other sediment types have been found: (1) peat, primarily on the south rim of the lake; (2) sand, in the nearshore region on the west side of the lake and (3) marl/rock, along an outcrop just north of the peat region (Fisher, Reddy, and James 2001).

Sediment cores were obtained from Lake Okeechobee in 1988, 1998 (Fisher, Reddy, and James 2001) and 2006 (Balance Environmental Management (BEM) Systems and University of Florida 2007; Yan and James 2012) at 170 sampling points (Figure 2). The thickness of mud sediments at the top of each core was measured. Fisher, Reddy, and James (2001) describe methods used to extract and analyse total nitrogen (TN), TP, iron (Fe), calcium (Ca) and carbon (TC) from the top 10 cm of each core (Table 1). These values along with their coordinate locations were used to create a number of geospatial models as described below.

Table 1. Summary table of nutrients (mg/g), mud thickness (cm) and elevation (ft).

		TC	CA	FE	TN	TP	Thickness	Elevation
1988	Min.	0.6	41	24	166	38.6	0	−1.27
	1st Qu.	17.75	7297	643	1541	220.7	0	0.31
	Median	107.8	54,377	1700	8621	664	0	2.6
	Mean	118.93	53,799	2695	9094	670.5	12.47	4.32
	3rd Qu.	151.8	85,154	4787	11,319	1076	18.88	8.14
	Max.	481.9	176,020	10,524	36,107	1708.1	66	15.44
	Skewness	1.46	0.48	0.61	1.38	0.16	1.26	0.66
	Kurtosis	1.89	−0.86	−0.56	1.68	−1.42	0.29	−0.91
1998	Min.	1.4	220	112	0	2.2	0	−1.27
	1st Qu.	25.93	10,396	1449	400	165.7	0	0.31
	Median	143.6	51,325	3506	8000	474.7	0	2.78
	Mean	145.06	60,705	7060	9363	650.2	11.17	4.39
	3rd Qu.	182.45	89,439	14,352	12,500	1149.8	13.75	8.15
	Max.	487.1	328,000	21,123	47,800	1793.3	74	15.44
	Skewness	1.06	1.47	0.60	1.47	0.30	1.49	0.61
	Kurtosis	0.35	2.46	−1.35	1.84	−1.48	0.99	−0.95
2006	Min.	1.15	114	39	70	16.64	0	−1.12
	1st Qu.	25.69	6813	716	610	141.39	0	0
	Median	140.77	64,349	2013	8170	397.42	1	0.99
	Mean	139.09	64,217	3429	8926	566.22	8.27	3.16
	3rd Qu.	175.81	94,162	6947	12,500	999.46	14.75	5.94
	Max.	490.55	347,790	8827	41,030	1285.31	51	14.58
	Skewness	1.10	1.68	0.38	1.15	0.13	1.28	1.06
	Kurtosis	0.69	3.83	−1.53	1.07	−1.71	0.65	−0.20

In addition, the lake bottom elevation was surveyed during July–September 2008. The 681,796 survey data points were integrated with a LiDAR (Light Detection and Ranging) survey data collected during 2007–2009 for the lakeshore areas. These data sets were combined to develop an accurate digital elevation model (500 feet by 500 feet) for the entire area of Lake Okeechobee (Figure 2). Elevation values for the 170 sample locations were extracted from the integrated elevation raster to be used in geospatial models.

GWR models

There are a number of assumptions associated with spatial statistical prediction models. Variables of interest must be stationary, which means that the statistical properties of the variables (e.g. mean, standard deviation and covariance) are the same over space and time for the area of interest. However, spatial data exhibit two properties: spatial dependence and spatial heterogeneity, and this makes it difficult to meet the assumptions and requirements of traditional (non-spatial) statistical methods. For example, geographical properties vary from place to place. A local spatial model can estimate temperature values at un-sampled locations, better than a global model, based on the values at nearby locations and known values of elevation (or other ancillary variables).

There are different strategies to deal with spatial autocorrelation in regression model residuals, such as re-sampling until the input variables no longer exhibit statistically significant spatial autocorrelation, or isolating the spatial and non-spatial components of each input variable using a spatial filtering regression method. For the traditional non-spatial OLS models, regional variation (non-stationarity) can be eliminated by defining/reducing the size of the study area so that the processes within it are all stationary. A better approach is to use methods that incorporate regional variation into the regression model such as GWR (Brunsdon, Fotheringham, and Charlton 2010; Fotheringham, Brunsdon, and Charlton 2002).

GWR is a spatial regression method developed by Brunsdon, Fotheringham, and Charlton (2010), Fotheringham, Brunsdon, and Charlton (2002), building on the works of Hastie and Tibshirani (1990) and Loader (1999). It can calibrate a multiple regression model that allows different relationships between the dependent and independent variables to exist at different locations in space (non-stationarity). GWR is a local spatial regression approach based on the ‘First Law of Geography’: everything is related to everything else, but closer things are more related to each other (Tobler 1970). As a local analysis technique, GWR tends to serve as a ‘microscope’ that amplifies details of the data that are otherwise hidden. GWR provides a local model of the predicted variable or process by fitting a regression equation to every feature in the data set. When used properly, these methods are powerful and reliable statistics for examining linear relationships. With these characteristics, GWR is recommended as an additional spatial interpolation technique that may increase interpolation accuracy. GWR assumes that relationships between the target variable and ancillary variables are not constant over space, and the error terms are independent and identically distributed (i.i.d.) across space. Therefore, the coefficients (β_ks) vary from location to location; the interpolation can then be expressed as

(1)

where y_i is the predicted value at the , β₀ is the model intercept, β_k is the coefficient of the kth explanatory variable, n is the the number of available ancillary variables, x_k = kth is the explanatory variable, k = 1, …, n, ε_i is the error associated with the ith observation, (u_i v_i) are the coordinates of the ith observation in space and is a realization of the function at point i.

It is assumed that each prediction is generated through a Gaussian or Gaussian-like spatial process in which the strength of the correlation among observations declines as distance increases. With such an assumption, it is possible to mimic the spatial process at any particular location in space, and generate a set of local predictions that are based on weighting of available observations according to the distance and the prediction location. The raster surfaces of model coefficient values created by GWR can be examined to evaluate regional variation in the model explanatory variables. The spatially consistent (stationary) relationships between the dependent variable and each explanatory variable and how they change across the study area can be evaluated as well as the coefficient distribution that shows where and how much variation is present.

GWR also can be used for prediction when it is applied to sampled data. First, a data set is specified with all of the explanatory variables for locations where the dependent variable is unknown. GWR then calibrates the regression equation using known dependent variable values from the input feature class, and creates a new output feature class with dependent variable estimates.

Weighting functions and bandwidth selection

The parameters estimated for a GWR model are dependent on the spatial weighting function and the bandwidth selected. The weighting function determines the importance of a local observation based on its proximity to the sample point (Figure 3). The bandwidth determines which nearby observations are considered when calibrating coefficients for a sample point. Bandwidths may be constant (fixed kernel) or variable (adaptive kernel) that vary by the location of a data point.

Figure 3. GWR weighting schemes and bandwidths. (a) Fixed weighting scheme. (b) Adaptive weighting schemes.

The weighted function for parameter estimation was developed by Fotheringham, Brunsdon, and Charlton (2002):

(2)

where the bold symbols are matrices, is an estimate of β and is an n by n matrix whose off-diagonal elements are zero and whose diagonal elements denote the geographic weighting of each of the n observed data for regression point i.

The weights are chosen such that those observations near the point in space where the parameter estimates are derived have more influence on the result than observations further away. Two schemes are used for the weight calculation: bi-square and Gaussian. In the case of the Gaussian scheme, the weight for the ith observation w_ij is

(3)

In the bi-square function

(4)

where d is the Euclidean distance between the location of observation i and location j, and h is the bandwidth.

With the fixed kernel-weighting scheme, the same shape and magnitude of the selected weighting functions are applied to every sample point over the space (Figure 3a). Fixed schemes might produce estimates with large variances where data are sparse, while masking subtle local variations where data are dense. Adaptive schemes adjust according to the density of data: shorter bandwidths are selected where data are dense and longer bandwidths selected where data are sparse. Adaptive kernels are often used when data are not evenly distributed (Figure 3b).

An optimal bandwidth (or nearest neighbours) can be selected by satisfying either least cross-validation (CV) score (the difference between the observed value and the GWR calibrated value), or least Akaike Information Criterion (AIC). Four GWR sub-models – Adaptive AIC_c, Adaptive CV, Fixed AIC_c and Fixed CV – are the different combinations of weighting functions and optimization criteria. In GWR, a corrected AIC (AIC_c) was used (Hurvich, Simonoff, and Tsai 1998), and it was first suggested for normal linear regression by Sugiura (1978). Hurvich, Simonoff, and Tsai (1998) demonstrated the small-sample superiority of AIC_c over AIC, and justified the use of AIC_c in the frameworks of non-linear regression models and autoregressive models. The AIC_c can be calculated using the following equation:

(5)

where n is the number of observations in the data set, is the estimate of the standard deviation of the residuals, S is the hat matrix and tr(S) is the trace of the hat matrix.

Cross-validation minimization can be used to search for an optimal bandwidth:

(6)

Model development

OLS models of observed data were first developed to find significant relationships between TP and other ancillary data. After one or more candidate regression models are identified using the ArcGIS 9.3.1 OLS statistical regression tool, the GWR model is applied for the selected regression models. The ArcGIS GWR statistical tool creates coefficient raster surfaces for the model intercept and each explanatory variable.

Corrected AIC_c was the measure chosen to determine goodness-of-fit. The AIC_c compares models of the same target variable and it contains a penalty for the complexity (degrees of freedom) of the model. The AIC_c provides a measure of the information distance between the fitted model and the unknown ‘true’ model. This distance is not an absolute measure but a relative measure known as the Kullback–Leibler information distance. AIC_c can compare the global OLS model with a local GWR model. Other model performance indicators were also calculated: the sum of the squared residuals (SSR), sigma (the estimated standard deviation for the residuals) and R².

Training and validation data sets were created with the Geostatistical Analyst extension of ArcGIS 9.3.1 by randomly selecting the 170 samples into a 75–25% split. ArcGIS OLS tool and Stepwise regression tools of R software were used to select the ‘best’ subset of predictors. Redundant (collinear) variables were removed using the variance inflation factor (VIF < 7.5). Two measures of model accuracy and precision were determined for validation points:

• Mean prediction error (ME):

(7)

• Root mean square prediction error (RMSE):

(8)

where l is the number of validation points, and are estimated and measured values at validation points , E is the expected value of the variable and is the variance.

In order to compare accuracy of prediction between variables of different types, the RMSE was normalized by the total variation. This is called the normalized root mean squared deviation or error (NRMSD or NRMSE):

(9)

where X_max and X_min are the maximum and minimum predicted values in the validation data points, respectively. Equation (9) determines the amount of global variation explained by the model. A value of NRMSE that is close to 40% means a fairly satisfactory accuracy of prediction (R-square = 85%). If NRMSE > 71%, then the model accounts for less than 50% of variability at the validation points (Hengl 2009). The validation results were calculated through a series of GIS processes using the ModelBuilder tool.

We also developed a histogram of errors at validation points and compared the errors estimated by the model and the true mapping error at validation points. This can help us detect outlier points where the errors are extreme compared to other locations. Scatter plots of predicted versus the measured values at validation points are used to derive a correlation coefficient to test the model’s ability. For comparison, OK, OLS and RK models were also examined and their errors were calculated.

In this study, the best GWR models were selected based on AIC_c, SSR (or sigma) and the adjusted R². If at least two of the SSR, AIC and adjusted R² values for Model A were better than those for Model B, while the remaining one (if any) was similar to that of Model B, then Model A was considered better. The model output feature class table includes fields for observed and predicted y-values, condition number (cond), Local R², residuals, explanatory variable coefficients and standard errors and can be examined for local errors and local collinearity.

Results

Mud thickness varied temporally and spatially. The maximum mud thickness for 1988, 1998 and 2006 were 66, 74 and 51 cm, respectively. The mean thickness declined slightly from 12.47 cm in 1988 to 11.17 cm in 1998, then dropped to 8.27 cm in 2006 (Table 1). The values of TP in sediments formed a bimodal distribution for each year sampled (Figure 4). TP varies in different zones, with the highest values in mud zones (>850 mg/kg), lower values in peat zones (250–450 mg/kg) and the lowest in sand and rocks zones (0–300 mg/kg) (Yan 2011). The 2006 TP values ranged from 17 to 999 mg/kg, with mean and median values of 566 and 397 mg/kg, respectively. The 1998 TP values ranged from near zero to over 1793 mg/kg, with mean and median values of 650 and 475 mg/kg, respectively. The 1988 TP ranged from 39 to 1708 mg/kg, with mean and median values of 670 and 664 mg/kg, respectively (Table 1). The TP concentrations of the three data sets show a second-order trend change along the north–south direction and relatively weak trend (first-order trend) in the east–west direction.

Figure 4. TP histograms (mg/kg) for 1988, 1998 and 2006.

Correlations (r) of the measured values show that TP is strongly and positively correlated with Fe (0.79–0.94), moderately and positively correlated with mud thickness (0.57–0.72), and weakly and negatively correlated with site elevation (−0.47 to −0.55) in all three data sets. Other parameters were not significantly correlated to TP (Table 2).

Table 2. Correlation coefficients and their significances of sediment nutrients (mg/kg), mud thickness (cm) and elevation (ft).

Year		Thickness	Elevation	TP	TN	CA	FE	TC
1988	Thickness	1
	Elevation	−0.5276	1
		<.0001
	TP	0.6046	−0.5472	1
		<.0001	<.0001
	TN	0.1255	0.1959	0.3703	1
		0.1222	0.0152	<.0001
	CA	0.4056	−0.6705	0.6460	0.1286	1
		<.0001	<.0001	<.0001	0.113
	FE	0.6540	−0.6864	0.7936	0.1460	0.6031	1
		<.0001	<.0001	<.0001	0.0718	<.0001
	TC	0.1583	0.1573	0.3382	0.9512	0.1187	0.1596	1
		0.0506	0.0522	<.0001	<.0001	0.144	0.0488
1998	Thickness	1
	Elevation	−0.5982	1
		<.0001
	TP	0.5715	−0.5474	1
		<.0001	<.0001
	TN	0.0301	0.2789	0.3612	1
		0.7131	0.0005	<.0001
	CA	0.2387	−0.5161	0.2611	−0.1659	1
		0.0031	<.0001	0.0012	0.0411
	FE	0.6932	−0.7038	0.8492	0.1469	0.2876	1
		<.0001	<.0001	<.0001	0.0709	0.0003
	TC	0.0675	0.1880	0.3533	0.9506	0.0286	0.1741	1
		0.4085	0.0204	<.0001	<.0001	0.7268	0.0319
2006	Thickness	1
	Elevation	−0.4678	1
		<.0001
	TP	0.7159	−0.4552	1
		<.0001	<.0001
	TN	0.1414	0.1452	0.37545	1
		0.0793	0.0714	<.0001
	CA	0.2918	−0.4579	0.3774	−0.0662	1
		0.0002	<.0001	<.0001	0.4134
	FE	0.7272	−0.4862	0.9374	0.3025	0.3603	1
		<.0001	<.0001	<.0001	0.0001	<.0001
	TC	0.1098	0.0925	0.3332	0.9616	0.0514	0.2556	1
		0.174	0.2522	<.0001	<.0001	0.5257	0.0013

Note: Upper values are the coefficients and bottom ones are the p values based on Pearson correlations.

Model selection and calibration

Two statistically significant OLS models with different independent variables were selected for TP concentration: total Fe alone, and mud thickness (Th) & elevation (Elev) together, respectively (Table 3). Each parameter was significantly related to TP independent of the others in the OLS models. The Global Moran’s I Index was used to calculate the spatial autocorrelation of the residuals.

Table 3. Statistically significant calibrated OLS models for Iron (Fe) and mud thickness and lakebed elevation (Th & Elev).

	1988		1998		2006
Predictors	Fe	Th & Elev	Fe	Th & Elev	Fe	Th & Elev
AIC	1580.78	1657.35	1563.97	1686.89	1461.83	1617.15
R^2	0.72	0.46	0.81	0.46	0.89	0.58
Adj R^2	0.72	0.45	0.81	0.45	0.89	0.57
F-Stat	287.33	47.38	489.51	47.69	907.33	76.11
K(BP)	2.01	3.00	1.56	1.51	0.03	3.45
K(BP)-Prob	0.16	0.22	0.21	0.47	0.86	0.18
JB	22.71	10.30	447.17	29.77	40.41	8.42
Sigma	246.07	342.80	228.58	390.23	146.04	287.40
Neighbours	15,795	15,795	13,665	13,665	14,729	14,729
Stationary	Yes	Yes	Yes	Yes	Yes	Yes
Residuals	Clustered	Clustered	Random	Random	Random	Clustered

For all three data sets total Fe as a predictor alone explained over 70% of the model variation while mud thickness explained over 45% of the variation (Table 3). The model residuals were either clustered (1988 and 2006 data sets) or randomly distributed (1998) as shown by the Jarque–Bera (JB) statistic. Koenker’s Breusch–Pagan (BP) value K is not significant, indicating the correlation between TP and total Fe is stationary in the study area. Correlation analysis shows that both mud thickness and elevation are statistically independent. The K (BP) statistic is not significant, indicating that the relationship between TP and thickness and elevation are also stationary in the study area.

GWR (TP vs. mud thickness and elevation) model

Four GWR sub-models (Adaptive AIC_c, Adaptive CV, Fixed AIC_c and Fixed CV) were developed using mud thickness and elevation as independent predictors. For all data sets the Fixed CV model had the best overall results (Table 4). For the 2006 data set, the Fixed CV model has the lowest errors (SSR, Sigma), better model performance (lower AIC_c value) and the model explained much of the variability (higher R²). This sub-model was further evaluated in the validation process. Its residuals were dispersed based on the Moran Index value (Yan 2011). For the 1998 TP data, all the GWR sub-models performed similarly with small differences and weakly dispersed residuals.

Table 4. Analyses of the GWR calibration model (TP vs. Th & Elev).

Year	GWR Models	Adaptive AICc	Adaptive CV	Fixed AICc	Fixed CV
2006	Neighbours	49	54	30,916	28,774
	Residual squares	4,738,304	4,895,445	4,211,677	4,010,590
	Effective number	17.73	16.03	25.51	28.09
	Sigma	221.85	223.54	218.17	216.068
	AICc	1569.80	1570.11	1571.12	1571.41
	R^2	0.78	0.77	0.81	0.82
	R^2 Adjusted	0.74	0.74	0.75	0.76
1998	Neighbours	66	66	44,943	42,884
	Residual squares	12,781,568	12,781,568	12,902,301	12,735,949
	Effective number	11.74	11.74	13.02	13.88
	Sigma	353.55	353.55	357.46	356.66
	AICc	1671.51	1671.51	1673.89	1673.93
	R^2	0.59	0.59	0.59	0.59
	R^2 adjusted	0.55	0.55	0.54	0.54
1988	Neighbours	32	27	30,777.37	21,021
	Residual squares	5,694,144	5,233,681	6,224,684	4,734,359
	Effective number	26.45	31.02	24.47	41.02
	Sigma	255.03	251.14	263.68	254.70
	AICc	1610.11	1611.91	1613.47	1625.79
	R^2	0.76	0.78	0.74	0.80
	R^2 adjusted	0.70	0.71	0.68	0.70

Note: In bold the most accurate models.

The model coefficient maps show very strong spatial gradients of the relationships to TP concentrations (Figure 5). In the 2006 data, mud thickness was positively correlated with TP, with a decreasing correlation from the eastern open water region of the lake to the western littoral region of the lake (Figure 5); site elevation was negatively correlated to TP in the central lake zone and positively correlated in most other regions of the lake. The 1988 and 1998 models had similar relationships (Figure 5).

Figure 5. Coefficient maps of the GWR model (TP vs. Th & Elev) calibration for 1988, 1998 and 2006.

GWR (TP vs. total Fe) model

For 1988 and 2006 data, all GWR sub-models performed better than the corresponding OLS model based on AIC_c and R². Among the GWR sub-models, CV models (Adaptive CV and Fixed CV) were better than the AIC_c models. The two CV models were very close in both AIC_c and R² values, but the Fixed CV has much lower SSR. The Fixed CV model was considered a better model for TP estimation (Table 5). For 1998 data, the Fixed AIC_c model performed the best. Total Fe was positively correlated with TP in all three data sets and the coefficients increased from east to west (Figure 6).

Table 5. Analyses of the GWR calibration model (TP vs. Fe).

	GWR models	Adaptive CV	Adaptive AIC	Fixed CV	Fixed AIC
2006	Neighbours	32	38	25,048	28,061
	Residual Squares	1,474,546	1,559,402	1,370,525	1,479,068
	Sigma	125.79	127.15	125.87	127.53
	AICc	1442.36	1442.05	1445.88	1445.18
	R^2	0.93	0.93	0.94	0.93
	R^2Adjusted	0.92	0.92	0.92	0.92
1998	Neighbours	114	46	2,583,313	38,985
	Residual Squares	5,721,282	4,772,090	5,851,677	4,754,339
	Sigma	228.39	219.36	228.58	217.70
	AICc	1566.66	1563.92	1566.19	1560.32
	R^2	0.82	0.85	0.81	0.85
	R^2 Adjusted	0.81	0.83	0.81	0.83
1988	Neighbours	13	38	21,434	24,337
	Residual Squares	2,256,157	4,669,867	3,403,058	3,809,536
	Sigma	189.93	220.30	205.77	210.14
	AICc	1582.82	1567.64	1564.55	1563.34
	R^2	0.91	0.81	0.86	0.84
	R^2 Adjusted	0.83	0.77	0.80	0.79

Note: In bold the most accurate models.

Figure 6. Coefficient maps of the GWR model (TP vs. Fe) calibration for 1988, 1998 and 2006.

GWR model validation and error analysis

Using elevation and mud thickness as independent variables, TP values were estimated on a 500 ft by 500 ft grid. The predicted values at the validation sites were compared to observed data. The residual errors for the 2006 validation data were normally distributed within a narrow region in the quantile-comparison plots (Yan 2011). Both residual errors of 1998 and 1988 were skewed with some large positive errors (underestimates) in the 1988 data. All scatter plots demonstrated moderate fits for the measured and estimated values with R² ranges between 0.67 and 0.69 (Figure 7).

Figure 7. Scatter plots of the GWR models (TP vs. Th & Elev) validation set for 1988, 1998 and 2006.

The coefficient values and intercepts change over the study area for Fe (Figure 5), and both values were used to calculate TP. All the residual errors were skewed distributions with higher positive errors (underestimates) for 2006 and higher negative errors (overestimates) for both 1988 and 1998 data sets (Yan 2011). Modelled versus predicted comparisons of the 2006 validation data were very good, with R² of 0.88 (Figure 8). With several extreme values removed from the 1988 and 1998 validation data, the comparisons of the modelled versus predicted were also good as well with R² of 0.87 and 0.86, respectively.

Figure 8. Scatter plots of the GWR models (TP vs. Fe) validation set for 1988, 1998 and 2006.

Model comparison and TP mass calculations using optimal models

In addition to GWR models, errors of OK, OLS and RK models were determined (Table 6). For all models evaluated, the GWR models (TP vs. Fe and TP vs. Th & Elev) had lower RMSE and the lowest ME for TP 2006 and TP 1988 data sets (Table 6). For 1998 TP data, OK has the lowest mean errors but higher RMSE than the OLS and GWR models.

Table 6. Comparison of TP model errors.

	Models	OK	OLS(Fe)	GWR(Fe)	OLS (Th, Elev)	GWR (Th, Elev)	RK
2006	ME	−27.38	31.19	−30.19	40.21	20.82	65.15
	RMSE	331.00	169.96	204.61	325.23	307.60	325.10
	NRMSE	26.09	13.40	16.13	25.64	24.25	25.63
1998	ME	−89.46	136.93	132.81	106.64	100.58	126.48
	RMSE	415.10	374.75	359.83	442.98	389.38	419.14
	NRMSE	23.18	20.92	20.09	24.73	21.74	23.40
1988	ME	−41.74	110.94	109.59	64.15	24.39	51.64
	RMSE	402.10	368.36	375.88	371.97	234.34	371.99
	NRMSE	0.24	0.22	0.23	0.22	0.14	0.22

Note: OK: Ordinary Kriging; OLS: Regression, Th: Thickness, Elev: elevation, RK: Regression-kriging; NRMSE: RMSE/(Max-Min); (In bold the most accurate models).

Based on model validation errors, the GWR (Fe) and GWR (Th and Elev) models were used for TP concentration prediction and mass calculation. For the GWR (TP vs. Fe) model, total Fe was estimated using OK for the whole lake bottom first. TP concentrations (a, c and e of Figure 9) and mass (b, d and f of Figure 9) were calculated using ArcGIS tools. For the GWR (Th and Elev) model, the mud thickness raster and the lake bottom elevation data raster were used to calculate the TP concentrations (a, c and e of Figure 10) and mass (b, d and f of Figure 10).

Figure 9. TP concentration (a, c and e) (mg/kg) and weights (b, d and f) (kg) estimated using the GWR model (TP vs. Fe) for 1988, 1998 and 2006.

Figure 10. TP concentration (a, c and e) (mg/kg) and weights (b, d and f) (kg) using the GWR model (TP vs. Th & Elev) for 1988, 1998 and 2006.

The results showed similar changes for TP mass: increase from 1988 to 1998 by 38–44%; decrease from 1998 to 2006 by 30–34%. From 1988 to 2006 (28-year period), there were small declines in TP mass: 1.0E + 06 kg or 2% reduction for OK and GWR (Fe) models, and 4.2E + 06 kg or 10% reduction for the GWR (Th and Elev) model (Table 7).

Table 7. TP mass (kg) and changes over time.

	OK model			GWR (Fe) model			GWR (Th + Elev) model
Year	Mass	Change	Change%	Mass	Change	Change%	Mass	Change	Change%
1988	4.25E + 07			4.20E + 07			4.43E + 07
1998	5.89E + 07	1.64E + 07	39%	6.04E + 07	1.84E + 07	44%	6.11E + 07	1.68E + 07	38%
2006	4.14E + 07	−1.75E + 07	−30%	4.10E + 07	−1.94E + 07	−32%	4.01E + 07	−2.10E + 07	−34%

Conclusions

GWR models provided the best predictions of TP in sediments for Lake Okeechobee because of the strong positive relationships between TP and total Fe, and moderate positive relationships between TP and mud thickness and site elevation. The small search windows used on local samples for estimation were an advantage of GWR models over global OLS models. OK models used only the weak spatial autocorrelation among TP data and led to the highest RMSE; their smaller mean errors are mainly due to smaller search windows used in variogram fitting. Regression-Kriging models use both weak spatial autocorrelation among TP values, weak correlation between TP and mud thickness and site elevation, and they performed better than OK models. These calibrated models may add new approaches to analyse other lakes or water bodies that are faced with similar challenges such as Lake George, Uganda and Lake Tai-Hou, China (Havens et al. 2007).

Mud sediments in Lake Okeechobee contain high amounts of P (Reddy, Sheng, and Jones 1995). Much of this mud has accumulated since the 1940s (Brezonik and Engstrom 1998). Brezonik and Engstrom (1998) found that these sediments consist of approximately equal amounts of organic, nonapatite inorganic-P and apatite inorganic-P. Thus, correlations of P to total organic carbon and calcium (found in apatite inorganic-P) in these mud sediments show the value of these ancillary measurements for developing robust models. Using the GWR technique, both organic carbon and calcium could be used to develop P predictions in Lake Okeechobee sediments. However, thickness of mud sediments, which contains an easily measured surrogate for organic content, and bottom elevation, which indicates where sediments will accumulate more readily (Jin, Ji, and Hamrick 2002), produce more accurate models than Ca and TOC.

Discussion

Sediments accumulate in the central deep regions of Lake Okeechobee (Brezonik and Engstrom 1998). This is due primarily to the wind-driven gyres that are set up on the lake and tend to deposit sediments in the deeper central areas of the Lake (Jin, Ji, and Hamrick 2002). Over time this material has collected in these regions, and has been mostly undisturbed as demonstrated from core dating using Pb²¹⁰ (Brezonik and Engstrom 1998) and the observance of sub-millimetre laminations in the mud sediments that disappeared only within the top few centimetres (Kirby, Hobbs, and Mehta 1994). The dramatic changes in the thickness and elevation coefficients for the GWR model (Figure 5) between the 1998 and 2006 period are indicative of the three major hurricanes – Frances, Jeanne and Wilma – that passed close to the Lake on 5 September 2004, 26 September 2004 and 24 October 2005, respectively. The hurricanes re-suspended a large amount of sediment and redistributed it throughout the lake (Jin et al. 2011).

The strong positive relationship between TP and Fe is related to the binding of inorganic phosphorus to oxidized iron in the upper few centimetres of sediment (Moore et al. 1993). Moore et al. (1993) suggested that the abundance of iron in the Lake Okeechobee system produces conditions that control P exchange by precipitation and absorption in the upper oxic layer of sediment rather than diffusion from the anoxic sediments below. These authors demonstrated the process by incubating sediment cores under alternating oxic/anoxic water column conditions. They found higher fluxes of P under anoxic conditions, suggesting that oxidized Fe was controlling the rate of this flux. This indicates that more soluble P could be sequestered under sediment layers with higher Fe. Because Fe is observed in regions with and without mud, it is a much better predictor than mud thickness. Fe also is better than elevation in predicting P because it is mobile and it affects the geochemistry of P.

Moore and Reddy (1994) used a similar framework to investigate the role of pH on P exchange from Lake Okeechobee sediments. They found that very acidic conditions (pH 5.5) resulted in SRP release. However, Lake Okeechobee pH is somewhat basic (approximately 8.2 – James, Smith, and Jones 1995) and the sediment pH is rather stable (7.0–7.5 – Moore and Reddy 1994). Thus, pH was not a useful parameter in these GWR models.

Between 1988 and 1998 the lake received an estimated net load of phosphorus of 2584 metric tons (Havens and James 2005); from 1999 to 2006 the estimated net load was 1949 metric tons (James and Zhang 2008). The changes predicted by the models are much greater and variable than these numbers derived from nutrient budgets. The tremendous variability among the sediment-sampled cores is also confounded by the small number of samples per area (less than 1 per square km). The large amount of error could account for the differences between the net loads to the lake and net change observed in the models. By sampling more intensively in areas where large changes in TP and ancillary data occur, model accuracy could be improved and thus improve estimates of change.

Sediment–water interactions are important to the health of Lake Okeechobee. Re-suspension of the mud sediments, especially during high lake water levels, can lead to reduced submerged aquatic vegetation, which negatively affects fish recruitment (Havens et al. 2005). Diffusive phosphorus fluxes, primarily from the mud sediments, produce an internal load that is equivalent to external loads (Fisher, Reddy, and James 2005). This can lead to delays of in-lake phosphorus in response to eternal loads (James and Pollman 2011). Given a program that measures sediment content of TP and ancillary data, accurate GIS models can be developed to evaluate spatial variability on these lakes. Such information could be used to manage these lakes through targeted chemical or dredging techniques (Søndergaard, Jensen, and Jeppesen 2003).

Acknowledgements

This research was supported by South Florida Water Management District. The authors wish to thank Susan Gray, Jun Han, David Unsell, Ming Chen, Kim O’Dell and Jenifer Barnes for their comments and suggestions to the previous version of this manuscript. The authors also wish to thank the three anonymous reviewers for their valuable comments and recommendations, which enhanced the manuscript significantly.