Eqn. 1
In this equation S indicates a smooth “thin plate” regression spline function, which is the MGCV packages default smooth function (Wood 2003). When the S term contains more than one variable, it is a multidimensional smooth that considers the interaction between the terms. Station is treated as a categorical variable. ε is a normally distributed error term. The fit of this function to our data indicated that across all stations, this variability with depth was statistically detectable, and that there was statistically significant station to station variability (R2 = 0.78 (overall model), p < 0.001 (for all interaction terms but one (pressure * station 3.3))). Particle abundance normalized to LISST size bins decreased as particle size increased (Fig. 3).
Particle Size to Abundance Relationship
At all stations, there was a negative power law relationship between particle size and particle abundance.  The slope of the power law distribution, which is the slope of the relationship between log transformed particle abundance and log transformed particle size, ranged at most stations and depths from -3.5 to -4. However, several depths at some stations had anomalously large negative particle size distribution slopes (Fig. 4). A general additive model of form shown in Eqn. 2 suggested that that there was statistically significant variability in the particle size distribution between depths, and that this relationship varied between stations (R2 = 0.167, p < 0.01).
           Eqn. 2
Total Particle Mass Patterns
Estimated total particle mass per liter of all particles > 1.2 μm ranged from 10 to 100 mg /L (Fig. 5). Calculated particulate matter concentrations were higher in the bottom sample than the surface sample at every station except 4.3 (OLS log(Mass) ~ Depth [Surface or Bottom, excludes Oxycline], F = 7.6, p = 0.02). At station 4.3 the sample taken in the oxycline had the highest particle concentrations, followed by the surface sample, and then the bottom sample. Particle concentrations estimated by LISST measurements were generally higher in the surface than in the bottom, except at stations 3.1 and 3.2. There was no detectable relationship between station latitude and observed particle mass (Ordinary Least Squares regression of form `log(Mass) ~ Latitude`; F = 0.001, p = 0.97).
Particle Mass to Size Relationship
Mass per particle increased with particle size, following a power law (Fig. 6).  The masses of particles of each size class were similar at each depth, ranging from about 10-9 mg/particle in the 1.2 µm class to about 10-3 mg/particle in the 500 µm class. There did not appear to be statistically significant differences between the slopes of the relationship between particle size and particle mass (Fig. 6). A linear model that explored interactions between size, station and depth on particle mass, of form `ln(Mass) ~ ln(Size) * Station * Depth`, where “ln” indicates the natural logarithm function, found that while there was a relationship between size and mass (p < 10-10), neither station, depth, nor any interaction between size, station and or depth had any statistically significant relationship to particle mass. However, a linear model of form `ln(Mass) ~ ln(Size) + Station` suggested that there was station to station variability in the intercept of the size to mass relationship (p < 0.01 for all stations, with the exception of stations 3.2 and 3.1 which had statistically identical intercepts). The `eemeans` package was used to compare the y intercepts of the size to mass relationship at each station. It was found that station 3.1 had statistically significantly lighter particles, adjusted for size, than stations 4.3 (difference = -2.5 +- 0.7 (1 standard error) log(mg/Particle), t-ratio = -3.86, p  = 0.012) or 5.1 (p = 2.5 +/- 0.7 log(mg/Particle), t-ratio = -3.79, p = 0.014). All other differences were found to be not statistically significant, after adjusting for multiple comparisons (FDR < 0.05). (Fig. 7).
Calculated Total Particle Mass Profiles
By combining the information from the mean particle size to mass relationship with the abundances of particles at each size, we were able to calculate expected particle mass throughout the water column at each station (Fig. 5; black circles). A general additive model of form shown in Eqn. 3 suggested that particle mass varied statistically significantly between depths (F = 9.4, p < 10-10), with all stations except 3.3 and 5.5 showing statistically significant deviations from the main profile (F >= 3.1 p <0.003 for all remaining stations).
                         Eqn. 3
Calculated total particle mass appeared to be related to, but was often an underestimate of, observed total particle mass (Fig. 5).
 

Discussion

Measurements of physical and chemical parameters (Fig. 1) showed depth profiles typical of previous measurements of the Chesapeake Bay at this time of year (Pritchard 1952; Murphy et al. 2011). The location of the stations arranged along the length of the Bay allowed for gradients to be observed. The salinity gradient in the Chesapeake is formed as colder, denser saline water enters the mouth of the Bay and flows northward, while warmer, less dense, freshwater enters from rivers and tributaries and moves south (Pritchard 1952). The density difference in these two layers forms a pycnocline, which was observed at all stations. The pycnocline blocks the vertical transfer of oxygen, creating the anoxic zones that were seen in most stations. Large anoxic and hypoxic zones form during summer in the Chesapeake Bay and were clearly seen in July when measurements were taken. Anoxic bottom waters have been shown to lead to increases in sulfide concentrations (Roden and Tuttle 1992), as seen in the sulfidic samples collected in stations 4.3 and 5.1.
As samples were collected at various times throughout the day, it is possible that changing tidal currents could have affected particle distribution, although tidal and wind current speed was not measured. The collection at station 4.3 occurred close to high tide, when stronger currents could potentially cause more particle resuspension. However, a study of tidal currents in the middle Chesapeake Bay found them too weak to resuspend bottom sediments significantly (Ward 1985).
Throughout the Chesapeake Bay, particle size distribution profiles displayed a power law relationship between size and abundance, with slope usually between -3.5 and -4, which is within the range of values seen in open ocean locations (Sheldon et al. 1972; Kostadinov et al. 2009; Cram et al. 2018). The slopes generally did not show much change with depth. This pattern is consistent with the findings of a study that size distribution does not change significantly with depth across the Atlantic Ocean (Gordon 1970), though it contrasts with measurements of an oxygen deficient zone in the Eastern Tropical North Pacific that found changes in the particle size slope with depth (Cram et al. 2022). The anomalous spikes of particularly negative slopes, seen especially in stations 4.3, 5.1, and 5.5, could indicate a lack of large particles in the oxycline, as the spikes occurred at approximately the same depth. These spikes could also be artifacts, perhaps induced by changes in salinity or temperature or introduced by the LISST’s electronics. Due to this possibility, we do not emphasize the spikes in our results and our overall conclusion remains that the particle size distribution slopes stayed mostly constant with changing depths.
While particle size spectra have been measured in the Chesapeake Bay, the particle size distribution slope is often not reported for the Bay (Schubel 1968; Schubel and Nelson 1973) or for estuaries in general. These measurements are more common in open ocean studies, thus providing values with which to compare our PSD measurements. However, it is expected that due to the much different conditions of open ocean water, there may be significant differences in PSD in estuaries that are yet unknown.  
The particle size to mass relationship also stayed consistent throughout stations and depths, with mass increasing and density decreasing in larger particles. This result suggests that particles of the same size may also have similar composition, leading to similar mass and density. The slopes of the size to mass relationship, or fractal dimension, at each station and sample depth (Figure 7) were similar to the values calculated in other particle studies. For instance, Fall et al. (2021) calculated fractal dimension in the York River as the size to density relationship with a bulk value of 2.25. Other studies have quantified fractal dimensions from the size to density or size to settling velocity relationship in the Chesapeake Bay (Sanford et al. 2004) and other estuaries and marine environments (Hill et al. 1998; Guidi et al. 2008; Jackson et al. 1997). These previous measurements of fractal dimension values for particles typically fall somewhere between 1.3 and 2.5, and the values for this study are on the low end of that range. Aggregation and disaggregation of particles affect their fractal dimension, with larger aggregates having lower fractal dimensions than small particles. Li and Logan (1995) found fractal dimension to decrease from 2.49 to 1.68 as particles coagulated during a phytoplankton bloom in a laboratory mesocosm. It is possible that collection methods in our study could lead to disaggregation of particles; however, the fractal dimensions’ consistency with other studies lend confidence to our observations.
Although particle size to mass relationships stayed consistent across stations and depth, total particle abundance and mass both varied by depth. Particle abundance profiles generally tracked total particle mass profiles (Fig. 2 and 5). The calculated total mass values from collected particles were consistently higher than the estimated mass based on LISST measurements, especially in station 4.3 (Fig. 5). This disparity could be caused by the assumption that the power law relationship between size and mass is the same at each station. Particle abundance and total particle mass both increased near the bottom of the water column, suggesting that the current is resuspending sediment from the floor of the Bay. In station 4.3 (Fig. 5), this increase at the bottom was seen in the mass from the samples (observation), but not in the LISST calculations, likely because the LISST measurements do not reach as deep as the sample, so the effect was not detected by the LISST. In the anoxic water below the pycnocline, particle abundance and total mass was lower than in surface waters. This scarcity of particle mass below the pycnocline suggests either fast removal of particles or low transport into this region from the more productive surface waters, as the pycnocline separates anoxic water from surface water. We argue that since particle remineralization is thought to be slow in anoxic water (Cram et al. 2018), it is likely the latter process, low in situ production and low flux from the surface that leads to the lower anoxic particle mass. Particle abundance and particle mass profiles diverged near the middle of the depth profiles, where abundance sharply decreased in most stations, but mass did not. The sharp decreases in particle abundance could be the same artifacts that may be seen in the particle distribution slopes. Again, these anomalies do not change the overall conclusion, with abundance and mass following similar patterns, highest near the surface and bottom of the water column.
This study was the first characterization of particle size to abundance and size to mass distribution across the length of the Chesapeake Bay as well as with depth. Overall, there was little variation by latitude, with particle size, abundance, and mass mostly following the same patterns at each station. This scarcity of variation pattern persists despite variations in particle source, with terrestrial sources in the north and in-situ production being more important in the south. We note that here we have only measured size and mass, and variability in other particle properties such as biology and chemistry may differentiate particles latitudinally. The sample stations covered the upper (CB3.1 and CB3.2) and middle Bay (CB3.3C, CB4.3, CB5.1) as well as the upper stream of the Rappannock Shoal (CB5.1) which falls just north of the lower Bay. The results suggest that within the areas studied, in the deep channel, it is factors other than latitude that lead to variability of particle characteristics. We note that the Estuarine turbidity maximum, characterized by higher particle loading and different dynamics than elsewhere (Malpezzi et al. 2013; Sanford et al., 2001; Keller et al. 2014), is north of our northernmost station CB 3.1 and therefore not included in this analysis. Significant differences were observed vertically, with particle mass and abundance higher just above the floor and low in the body of the anoxic layer. This low total particle mass in the anoxic water suggests that the influx of particles from the photic zone sinking and mixing into the anoxic water occurs more slowly than carbon removal from the anoxic layer, which would occur by particles settling into the sediment or by decay and remineralization. The low particle abundances could thus indicate the presence of particle remineralization despite the lack of oxygen, because if remineralization was absent, we would expect carbon accumulation in this region even if particle flux into the region was low. A low input of carbon into this region would also suggest carbon limitation of microorganisms in the anoxic region.

Conclusion

This analysis of particulate organic matter provides data for particle size distribution and particle mass at surface and bottom depths across various stations in the mesohaline region of the Chesapeake Bay. This study was the first to analyze particle distribution at multiple locations in the Bay. Generally, particle size/abundance and size/mass relationships were similar between stations and depths. Particle abundance and mass mostly followed similar patterns to each other, decreasing in the anoxic zones, with an increase near the bottom of the Bay. The results show the influence of depth on particle distribution, while patterns stayed consistent throughout station latitudes at the time of sample collection in July.