Re-Examining Field-Surveyed Variations in Elevation and Soil Properties with a 1-m Resolution LiDAR-Generated DEM ()
1. Introduction
The purpose of this study is to re-examine 1977-dated soil properties for a hummocky farm field near Hartland in New Brunswick (NB) as published by [1] . In this study, the elevations across the field were surveyed along a 5-m grid. Soil textural, morphological, and chemical properties were determined along a 25-m grid, and also along 5 × 5 m grid each at an uphill and a downhill location. The properties so analyzed refer to plough layer depth (Ap), soil texture, CF %, total soil organic carbon (TOC, i.e., particulate combined with non-particulate organic matter), soil moisture content (MC), water retention at saturation point (SP), field capacity (FC) and permanent wilting point (PWP), pH, and soil extractable Ca, Mg, K, S, P, NO3-N, NH4-N, Fe, Mn, Cu, Zn, and Cs137. This re-examination was facilitated through:
1) the availability of the 1-m resolution LiDAR DEM generated in 2017 [2] , and
2) the metric cartographic depth-to-water delineation index DTW formulated in [3] .
This index can be used approximate the areal extent of very poor to poor, imperfect, moderate, well and excessive soil drainage next to permanent open water bodies and channels as DTW varies from 0 to 10, 25, 50, 100 cm and more, respectively [3] . For forested areas, permanent flow channels generally require at least 4 ha of upslope flow-accumulation areas (FA) for end-of-summer surface flow [4] . Seasonally affected flow effects on channel-adjacent soil properties can also be analyzed by reducing the upslope FA requirement to, e.g., 1, 0.25 and 0.1 ha. The purpose for doing so refers to demonstrating how local uphill to downhill water flow and retention patterns affect measured within-field soil property variations.
2. Methods
2.1 Study Area
The field described by [1] is 175 × 425 m in width and length, and is approximately located 1 km one km east of Harland, NB (Figure 1; 46˚18'23''N 67˚30'41''W).
The geological surface deposit which is <2 m deep at this location refers to a sediment-derived loamy lodgment till with minor calcareous content. The soil that developed in this till is classified as an Orthic Humo-Ferric Podzol within the Carleton soil association. Cultivation - started 120 years ago - transformed the original mound-and-pitted forest soil underneath northern tolerant hardwoods to a smoothed surface with interrupted soil layer sequences. Across the surveyed field and beyond, intensified crop management since 1950 including potato cropping induced slope-dependent soil erosion coupled with soil re-deposition in depressions, at a rate of 22 to 53 tons/ha/year. Mean annual precipitation amounts to 1096 mm, with 796 mm from May to September. The mean monthly May to September temperature is 14.9˚C. The mean annual air temperature is 4.0˚C.
![]()
Figure 1. Locator map for the 1977 field survey in Harland, New Brunswick (46˚18'23''N, 67˚30'41''W).
2.2. Soil Analysis
The surface soil (primarily A-layer) was sampled (augured) and analyzed along the 25 m and 5 × 5 m grid for soil texture (hydrometer method without OM removal), CF content, and TOC concentration (combustion method using a Leco CNS-1000 analyzer). Also determined were the soil pH (1:1 water (H2O)) and Mehlich-3 extractable Ca, Mg, K, S, P, Fe, Mn, Cu, and Zn [5] . The Mehlich-3 extract formulation is composed of 0.2 mole/liter (M) acetic acid (CH3COOH), 0.25 M ammonium nitrate (NH4NO3), 0.015 M ammonium fluoride (NH4F), 0.013 M nitric acid (HNO3), and 0.001 M ethylenediaminetetraacetic acid (EDTA). Additional determinations involved:
1) 0 to 15 cm depth soil moisture levels on 11 July (MC1) and 27 August (MC2) 1997 using time domain reflectometry.
2) Calcium chloride (CaCl2) extractable NH4-N and NO3-N.
3) Cs137, using a Tennelec germanium crystal gamma radiation counter.
4) Soil saturation point (SP), field capacity (FC) and permanent wilting point (PWP).
All these determinations were done at and obtained from the Potato Research Centre of Agriculture and Agri-Food Canada in Fredericton, NB.
2.3. GIS Analysis
The GIS analysis was conducted with ArcMap using the 2017 LiDAR-generated 1-m DEM and Tarboton’s D8 algorithm. Doing so generated the slope, filled, flow direction, flow accumulation, and flow channel rasters [6] . The latter were classified into flow channel networks with >4, 1, 0.5, 0.25, and 0.1 ha upslope flow accumulations for flow initiation. The slope and reclassified flow-channel raster were used to determine the cartographic cost-distance derived DTW (in m) so that the >4, >1, >0.5, >0.25, >0.1 ha DTW classifications would respectively represent DTW at the end of summer (>4 ha), following storm and rainfall events in summer (>1, >0.5, >0.25 ha), and at intense snowmelt times (>0.1 ha, Figure 2). Also determined was the Topographic Position Index (TPI, [7] [8] ) using mean 25-m annulus elevations as TPI = 0 m reference. The 1977 field-surveyed point data were subsequently supplemented with their corresponding 2017 DEM, DTW and TPI extracted values using the Multipoint Extraction tool.
2.4. Statistical Analysis
The combined and GIS generated point data were summarized using basic statistics followed by correlation and simple and multivariate linear regression analyses. The resulting correlation matrix was factor analyzed to reveal how the variables so summarized either relate or differ from each other by way of an explanatory three-factor pattern. The regression analyses were done to determine 1) how the field-surveyed and LiDAR-generated elevations vary in detail, and 2) how Sand %, Silt %, Clay %, CF %, SOM %, NO3-N, NH4-N, Soil Moisture Content, Field Capacity and Permanent Wilting Point are affected by topographic position in the form of the DTW variables as defined by their minimum upslope flow accumulation areas.
3. Results and Discussion
3.1. Data Summary
The data associated with the field-surveyed variables are presented in Table 1, together with their units, averages, standard deviations, and minimum and maximum values. The correlations for most of these variables are listed in Table 2.
3.2. Field-Surveyed (1977) versus LiDAR-Registered (2017) Elevations
The summary in Table 1 reveals that the elevation changes from 1977 and 2017 were generally minor after accounting for the 84 m difference in resetting the 1977 to the 2017 zero-elevation reference level. Using the 2017 DEM as the only predictor variable for the 1977 elevation data produced a regression coefficient of <1 (i.e., 0.941, Table 3, Analysis A). This means that the survey plot was somewhat flatter in 1977 than in 2017. Adding DTW > 4 ha and then Silt % to the regression analysis (Table 3, Analyses B and C):
![]()
Figure 2. Survey grid and cartographic DTW associated from A to D, with the LiDAR-DEM derived flow channels with >4, 1, 0.25, and 0.1 ha upslope flow-accumulation areas, overlaid on the hillshaded DEM, respectively. DTW grades from <10 cm (dark blue) to 1 m (light blue) deep.
![]()
Table 1. Statistical summary of the field-surveyed variables with unit, mean, standard deviation, and maximum and minimum values.
DTW, Slope, TPI derived from 2017 DEM; Slope: focal 5 m circle mean; TPI = 2017 DEM - mean 25 m 2017 DEM annulus.
![]()
Table 2. Correlation matrix for most of the variables listed in Table 1. Significant regression coefficients (<−0.300 or >0.300) are highlighted in gray.
1) produced negative coefficients for these variables,
2) increased the R2 values to 0.993 and 0.995,
3) rendered the 2017 DEM coefficient to become 1.010 and 1.002,
4) reduced the RMSE value of the residuals from 0.249 to 0.190 m.
The numerically flatter elevation profile across the 1977-surveyed field is in part related to 2017-DTW and 1977-Silt % adjustments. In detail, the 2017 DTW adjustments render the 1997 hilltop locations some smoother while the Silt % adjustments render the 1977 downhill locations somewhat deeper. Further adjustments towards a residual RMSE of 0.162 m were obtained by regressing the resulting 1977 DEM residuals (Table 3, Analysis D) against the 2017 Slope and 2017 TPI variables. The positive 2017 Slope coefficient suggests that some of the 1977 elevations were slightly higher along the steeper slopes than in 2017. The negative 2017 TPI coefficient implies that some of the knoll elevations were slightly less pronounced in 1977, and some of the depressed areas - where TPI < 0 m - were slightly more filled in 2017 than in 1977. Figure 3 illustrates the extent to which some of the Analysis C residuals correspond to the underlying 2017 DEM derived Slope % and TPI rasters.
![]()
Table 3. Multivariate regressions(A, B, C) with the 2017 LiDAR DEM, DTW > 4 ha and 1977 Silt %as 1997 DEM predictor variables, followed by analyzing the resulting 1977 DEM residuals using the 2017 determined Slope and TPI variables (Table 1) as independent variables.
![]()
Figure 3. Analysis C (Table 3) residuals overlaid on Slope % (left) and TPI (right), together with the corresponding actual versus best-fitted scatterplot (center) and the 2017 DEM - derived flow channels with >0.25 ha upslope flow accumulation (white lines).
3.3. Regression Analyses: Sand, Silt, Clay and Coarse Fragment
Plotting Sand, Silt, Clay, and CF % versus DTW > 4 ha showed that Sand % increased but Silt and Clay % decreased with increasing DTW > 4 ha (Figure 4, left). The decreasing Clay % content with increasing DTW > 4 ha relates directly to the upland Si % loss, i.e., clay displacement did not occur across the field. Hence, the higher lying areas were found to be coarser and sandier than the lower less well-drained areas. This would likely be due to natural and recurring cropping-induced upland-to-lowland silt-displacing soil erosion. In this regard, Figure 4 (right) shows how surveyed CF % follows the underlying DTW > 4 ha pattern. This is further illustrated in Figure 5 where the dotted 1977 CF % pattern is more closely aligned with the 2017 DEM derived DTW > 4 ha pattern than with the 2017 DEM pattern.
The very coarse (vc), coarse (c), and medium (m) sized fine (f) and very fine (vf) fractions of Sand also increased with increasing DTW> 4ha, but with the trend decreasing towards finer grain size such that vcSand > cSand > mSand and no discernable DTW > 4ha for vfSand (Figure 6, left). Testing to which extent the DEM-generated patterns for DTW > 4, > 1, > 0.25, and > 0.1 ha were related to the field-determined Silt % led to the regression results in Figure 6 (right). The corresponding R2 values change in the order:
DTW > 0.1 ha, R2 = 0.286; DTW > 0.25 ha R2 = 0.255; DTW > 1 ha, R2 = 0.329; DTW > 4 ha, R2 = 0.372.
This means that the field-assessed Silt % variations are best expressed by the slope-affected cost distance between each survey point and its closest > 4 ha down-stream location.
![]()
Figure 4. Sand, Silt, and Clay % (left) and CF % (right) versus DTW > 4 ha, all with regression equations.
![]()
Figure 5. Surveyed CF % overlaid on the hillshaded 2017 DEM grid (left) and cartographic 2017 DTW > 4 ha grid (right). Also shown: 2017 DEM-derived flow channels with >0.1 ha upslope flow accumulation areas.
![]()
Figure 6. Left: decreasing trend from very coarse (vc) to coarse (c) and fine medium (m) Sand % fractions versus DTW > 4 ha upslope flow accumulation. Right: Silt % versus DTW (m) along the 2017 DEM derived flow channels with DTW > 4 ha, >1 ha, >0.25 ha and >0.1 ha upslope. All with regression equations.
3.4. Regression Analyses: Soil Moisture Content and Retention
Field capacity (FC %) and permanent wilting point (PWP %) increased significantly with decreasing DTW, with best results obtained using DTW > 0.25 ha as independent variable (Figure 7, left). In contrast, the soil saturation point was not so affected, i.e. SP % = 50.8 - 0.039 DTW > 0.25 ha (in m); R2 = 0.001. In general, SP varies with soil bulk density, while FC and PWP increase with increasing soil organic matter and as soil textures become finer [9] . This suggests
![]()
Figure 7. Left: Field Capacity (FC %), Soil Moisture (MC1 % and MC2 %), and Permanent Wilting Point (PWP %) versus DTW > 0.25 ha, in m. Right: NO3-N and NH4-N versus log10 (DTW > 0.25 ha, m).
that the bulk density of the soil was not affected by topographic position, but FC and PWP would have been influenced by increased Silt % and Clay % at the lower DTW location. The soil moisture determinations MC1 and MC2 determinations for 11 July 1977 and 27 August 1977 were - in terms of dryness - closer to PWP than to FC, with some of the MC1 determination trending higher at low DTW levels. This was not the case for MC2 on 27 August 1977.
3.5. Regression Analyses: NO3-N versus NH4-N Retention
Figure 7 (right) reveals a significant trend towards increasing NO3-N levels with decreasing DTW > 0.25 ha. This relationship becomes even more significant by noting that four of the five NO3-N > 12 mg/g levels occurred close to the tracked-DTW > 0.25 ha flow paths in Figure 2. Adjusting the log10 (DTW > 0.25 ha) values for these points to 1 cm and deleting the remaining outlier modified the best-fitted regression result in Figure 7 (right) to become NO3-N = 5.12 - 2.5 log10 (DTW > 0.25 ha); R2 = 0.408. In contrast, the corresponding NH4-N pattern remained low with no significant DTW trend.
3.6. Influences of Topography and Other Factors on the Surveyed Soil Properties
Factor analyzing the correlation matrix in Table 2 revealed three factors that account for 46.9% of the total correlation variance. The resulting polygonised factor-to-factor association pattern is presented in Figure 8, showing Factor 2 versus Factor 1 and Factor 3 versus Factor 2. The Factor that is not represented along the x- and y-axes appears as the polygon at or near the center. In terms of the total variance represented by the correlation matrix, F1 accounts for 26.0%, and can be interpreted as a Soil Loss Factor with its positive loadings for DTW, CF % and Sand % content, and its negative loadings for moisture (MC1 %) and Silt %. Factor 2 accounts for 11.6 % of the total variance and can be interpreted as a Soil Cropping factor, with its Mehlich-3 > 0.5 loadings referring to Ca, Mg, K, P, S, NO3-N, and Mn. Factor 3 accounts for 9.3% of the total variance and can be interpreted as a organo-metal complexation (SOM, or TOC) factor due to the polygon-represented loadings involving TOC, Ap, Fe, Zn, Cu, and Cs137 [10] .
Note that there is an overall upland-to-downhill drift of the positive and negative F1 loadings for F2. This is likely due to persistent uphill-to-downhill transfer of water and soil sediments. Also note the positive Mn loading on Factor 2. This could be due to Mn applications intended to control common scab proliferations (Streptomyces scabies; [11] ). The positive TOC-linked F2 entry for Cu could be due to foliar Cu applications intended to control occurrences of potato blight (Phytophtora infestans), scab and black-scurf inducing Rhizoctania, and potato-damaging nematodes ( [12] [13] ). Similarly, the TOC complexed Zn loading to F3 could be due to Zn applications intended to improve tuber yields [14] . The negative < - 0.25 F3 loadings for NO3-N and NH4-N in the F3 versus F2 plot of Figure 8 reflect the tendency of soil organic matter to retain ions in the following order:
NO3-N << NH4-N < K ≈ Na ≈ Cs137 < Mg ≈ Ca < Mn ≈ Zn ≈ Cu < Fe.
3.7. Multivariate Analysis
The multivariate analysis results in Table 4 serve to elaborate on the Factor 1, 2, 3 patterns for the chemical soil properties in Table 1 in quantitative terms, as follows.
![]()
Figure 8. Factor analysis: Soil Cropping Factor (F2) versus Soil Loss factor (F1), left; Soil Organic Matter Factor (F3) versus Soil Cropping Factor (F2), right.
![]()
Table 4. Multivariate regression results for the soil chemical variables listed in Table 1.
1) Mehlich-3 extracted Mg is highly correlated with Ca (Equation (1)), possibly due to the presence of calcareous soil parent materials (e.g., the Carleton Forest Soil Association) and/or dolomitic Ca/Mg carbonate applications [15] .
2) Mehlich-3 extracted Ca increases with soil pH and P but decreases with increasing DTW > 4 ha (Equation (7)). Part of this would be due to pH-elevating Ca carbonate applications. Also, the flow of water-soluble Ca would enrich extractable Ca and Mg at low DTW > 4 ha field locations.
3) The elevating Ca effect on pH effect can also be noted with Equations (4) and (9). In general, increasing the soil pH:
a) facilitates increases in P availability;
b) compensates for soil-acidifying Ca, Mg, K, and NH4-N root uptake;
c) reduces the possibility of low-pH induced P fixation, and root-damaging Al (aluminum) and Mn mobilizations [16] [17] .
4) Equations (5) and (9) reflect that S applications not only involve pH-neutral CaSO4 (gypsum) and K2SO4 for adjusting S deficiencies, but also elemental S applications to enforce downward scab-eliminating pH adjustments [18] [19] [20] .
5) Equations (8) and (10) suggest that increasing NO3-N and NH4-N levels would in part be due to NO3-N and NH4-N applications, possibly involving KNO3, NO3NH4 and/or urea. Figure 7 (right) indicates that NH4-N content << NO3-N content. This suggests that applying NH4-N or urea as N fertilizer was likely not practiced due to NH4-induced a) soil acidifying nitrification [21] , b) greater potatoes tolerance for NO3-N than for NH4-N [22] , and c) NH4-N induced K displacement from soil cation exchange sites (Equation (10)). The slight NH4-N increase with increasing DTW > 4 ha (Equation (10)) corresponds with lower denitrification rates on well-aerated upland field locations [23] . The more significant NO3-N increase with decreasing DTW > 4 ha (Figure 7, right) would be due to uphill-to-downhill NO3-N leaching from coarse-textured soils with low anion retention capacities [24] .
6) As per Equations (8), (9), and (14), K increases with increasing NO3-N but decreases with increasing NH4-N. The former likely relates to KNO3 and K phosphate applications, while the latter would be due to K-induced displacement of exchangeable NH4-N. Altogether, Mehlich-3 extracted K is affected by four variables, namely pH, NO3-N, P, and DTW > 4 ha (Equation (14)).
7) Mehlich-3 extracted P would not only increase by way of Ca phosphate applications [25] , but is also seen to increase with increasing Mn and Sand % content (Equation (4)). The significant contributions of Mn to P (Equation (6)) could be due to Mn phosphate applications. The decrease in extractable P with decreasing Sand % could be due to lower P-fixing Fe content along the lower areas of the field [26] [27] .
8) Mehlich-3 extracted Cu, Fe, Mn, Zn, Cs137 are linked to one another via Factor 3, but their quantitative dependencies are element specific. In detail:
a) Mehlich-3 extracted Zn is primarily related to Mehlich-3 extractable Cu but weakly so with increasing DTW > 4 ha (Equation (2)).
b) Mehlich-3 extracted Mn and Cu are in part quantified by Mehlich-3 extractable Fe (Equations (6) and (11)), with Cu also increasing with increasing Ca, P, and Sand % content.
c) Mehlich-3 extracted Cs137 increases with Mehlich-3 extracted Cu, soil C, plough player depth and Silt % (Equation (12)). Hence, soil C and Mehlich-3 extracted Cs137 increase slightly from the upland to the lowland field locations.
The Mehlich-3 extracted Cu, Zn and Cs137 fractions are also indirectly related to increasing DTW 4 > ha in a positive or negative sense, i.e., positive for Zn (+) via the positive relationship between Zn and DTW (Equation (2)), negative for Cs137 via decreasing Silt % with increasing DTW (Figure 4), and negative for Cu via increasing Sand % with increasing DTW (Figure 4).
9) Mehlich-3 extracted Fe increased with increasing Mehlich-3 extracted Cu and Mn but decreased with increasing NO3-N and pH (Equation (13); [28] . In general, increasing pH leads to decreased Fe hydroxide solubility. Decreased Fe extractability with increasing NO3-N in low-lying and less aerobic field locations would therefore be due to loss of redox-solubilized Fe.
10) While Cs137 is associated with TOC (+), K (+), Cu (+), and Silt % (+) according to Equation (12) in Table 4 (R2 = 0.546), it is also higher at low DTW locations, i.e.,
CS137 = 1641 - 291 log10 (DTW > 0.25 ha); R² = 0.111.
In contrast, TOC does not vary significantly with DTW, i.e.,
TOC = 2.23 - 0.009 (DTW > 1 ha); R2 = 0.010.
The CS137 vs. DTW trend therefore contrasts the lack of a TOC vs. DTW trend. This difference is likely due to observations that Cs137 binds more strongly to mineral surfaces than to organic matter [29] . Hence, the downhill CS137 trend is likely due to erosion-induced uphill-downhill silt transfer Figure 9, [30] .
11) According to Table 2 correlations, TOC increased with increasing Silt, Clay, NO3-N and NH4-N content but decreased with increasing Sand and CF content. Among these, the increase with Silt % and the decrease with NO3-N are more significantly related to TOC (Equation (3)). Generally, TOC increases downhill in silt-accumulating depressions across hummocky fields [30] [31] [32] [33] , but - in this survey - uphill to downhill TOC did not vary significantly (Figure 10), i.e.,
TOC % = 2.20 − 0.018 log10 (DTW > 0.25 ha); R2 = 0.002 and
TOC % = 2.00 + 0.0039 DEM; R2 = 0.0018.
In part, this could be due that TOC % remains low due to increased downhill Silt % and Clay %. To that effect, TOC only varies from 1.4% to 3% while Silt and Clay increase from uphill to downhill by about 6 and 3%, respectively (Figure 4).
12) The effect of topography on in-field soil property variations was further detailed by [34] based on 774 survey points spread across New Brunswick. That study summarized the trends so observed in terms of three topographically related Factors, as follows:
a) The dominant Factor referred to the increasing Sand (+), CF (+) and Soil Organic Carbon (SOC, +) in association with decreasing Silt (−), Labile N (−) and Particulate Organic Carbon (POC, −). These trends are similar to the above results, except that SOC (≈4 POC) should be increasing with increasing Silt %, as is the case for this study (Table 2, Equation (3), Figure 10).
b) The second Factor refers to the uphill decreasing effects on pH (−) and labile N (−), as is the case in Table 2 with respect to pH and NO3-N.
c) The third Factor refers to an increasing effect of curvature on soil resistance to penetration (+), and CF % (+). Based on curvature across eroding knolls, CF % (+) can be expected to vary with Sand % (+), Silt % (−) and DTW (+). Adding curvature as independent variable to the analysis of the 1977 survey results did not affect the results reported above.
![]()
Figure 9. Scatterplots of Cs137 versus (a): log10 (DTW > 0.25ha, m); (b): Cu; (c): Sand; (d): Silt; (e): Clay; (f): TOC.
![]()
Figure 10. Scatterplots of TOC % versus (a): Silt %; (b): Sand %; (c): NO3-N; (d): P; (e): CF %; and (f): log10 (DTW > 0.25 ha).
4. Conclusions
In summary, it is now possible to quantify erosion- and water-flow induced changes in elevations and surveyed soil properties across fields, at high-resolution. In detail, the regression results in Table 3 suggest that the point-generated elevations were somewhat smoother in 1997 than in 2017. This change was likely due to continuing soil erosion, which would further expose upslope rocks and other coarse fragments and lowering up-field soil organic matter while deepening silt-enriched flow channels in the downslope locations.
Factor analyzing the variables in Table 1 revealed three variation-controlling factors. Factor 1 links the uphill-to-downhill pattern for sand, silt, clay, and CF to the cartographic DTW index. Factor 2 refers to periodic N, Ca, Mg, K, S, Mn and P soil amendments, while Factor 3 reflects the retention of heavy metals such Fe, Zn, Cu and Cs137 by SOM (or TOC), and includes increased NO3-N levels in low-lying areas due to low SOM anion retention.
The variables that are directly or indirectly affected by DTW refer to:
1) Sand %, Silt %. Clay % and CF %; these increased with increasing DTW > 4 ha;
2) MC1, FC, PWP; these increased with decreasing DTW > 0.25 ha;
3) NO3-N, Ca, Mg, Cu, Cs237; these decreased with increasing log10 (DTW > 0.25 ha) pattern.
The cumulative effects of soil erosion and retention manifest themselves at the DTW > 4 ha scale. In contrast, the cumulative effects of water flow and retention manifest themselves at the DTW > 0.25 ha and log10 (DTW > 0.25 ha) scales, with the former and latter pertaining to soil moisture content and Mehlich-3 element concentrations, respectively. In contrast, TOC % was not related to DTW in this study, but the reported uphill-to-downhill TOC % variations [30] [32] [33] likely vary at the DTW > 4 ha scale as well.
In the absence of locally repeated field surveys, one-time soil property assessments such as above do not lend themselves for quantifying cause-and-effect relationships related to sequential crop management actions. Nevertheless, the analytical results and related correlations so derived are consistent with 1) quantifying topo-induced soil property changes and 2) potato-cropping recommendations with respect to regular and/or intermittent N, P, K, Ca, Mg, S, Fe, Mn, B, Cu, and Zn applications as summarized (see, e.g., [35] ).
Acknowledgments
The preceding work is based on the Agriculture and Agri-food Canada field survey data centered on the Hartland potato field, with many thanks to Lien Chow for making these data available. Also much appreciated is the free GeoNB registered access to the 1-m resolution LiDAR DEM for New Brunswick, and financial assistance received from the New Brunswick Agriculture Department for crop suitability mapping.