Abstract

Background and Aims. Previous work in Australia has demonstrated the value of data-driven approaches to terroir analysis but, like other terroir research, focussed predominantly on the natural resources (soils, topography, and climate) on which winegrowing depends. In only very few cases have metrics of production performance also been considered. In this study, focussed on the Marlborough region of New Zealand, we integrated data pertaining to vineyard performance with biophysical data (soils and climate) describing the conditions under which grapes are grown to give a more holistic indication of regional-scale variation in the terroir of the Marlborough production system. Methods and Results. Digital map layers describing variation in climate, soil properties, and the yield and harvest date of Sauvignon Blanc (Vitis vinifera L.) were assembled and analysed for similarity in their patterns of spatial variation over six vintages (2014–2019) using k-means clustering. The results suggest that the Marlborough region has a characteristically variable Sauvignon Blanc production with crop phenology and harvest date strongly influenced by variation in temperature, and yield variation impacted by soil properties. Spatial variation in seasonal rainfall did not appear to impact on vineyard performance. Importantly, the Wairau and Awatere valleys which, hitherto, have been considered together as parts of a single Marlborough region, are shown to be distinct. Conclusions. This analysis is strongly suggestive of the Marlborough terroir being variable at the within-region scale. It also lends weight to the idea that estimates of vineyard performance in some parts of the region may be used to predict performance in others. Significance of the Study. The results have potentially important implications for the management of both vineyard operations and winery logistics, for wine marketing and for whole-of-industry planning around expansion or contraction. The methods used are free of any bias introduced to many previous studies of terroir zoning through adherence to historical or geopolitical boundaries, expert opinion of wines, and other heuristics.

1. Introduction

The development of the Marlborough grape-growing region, located in the northeast corner of the South Island of New Zealand, has occurred over the past 50 years. Starting from initial plantings in the early 1970s, largely in the central Wairau Plains, the region expanded to first cover the southern valleys of the Wairau, before progressively increasing in the mid-2000s into the Awatere Valley to the south, towards the coast of the Wairau Plains from the early 2010s, and then more recently to the upper valleys of the Wairau and Awatere valleys. The current area planted is approximately 28,883 ha (71% of the total New Zealand vineyard area), with 23,290 ha (or 81%) planted to Sauvignon Blanc (Vitis vinifera L.) predominantly as a single mass-selected clone of UCD1 (https://www.nzwine.com/media/21915/1-vineyard-report-2022.pdf). Marlborough Sauvignon Blanc was first recognised as a distinctive style, reflecting the terroir of the region, in the early 1990s when success, particularly in the United Kingdom market, resulted in increasing consumer demand; it is now a recognisable international wine style.

Terroir is a multivariate concept which describes the interaction between the conditions under which grapes are grown and wine is made and the sensory and chemical attributes of wine [16]. Thus, whilst there is no literal English translation of terroir, it can be considered to encompass many factors—soil, topography, climate, landscape, biodiversity, and management [7]—which give a wine its “sense of place” [8]. In a New Zealand context, terroir aligns closely to the Māori concept of tūrangawaewae [911]—literally, ‘a place to stand.’

Whilst management of both grapegrowing and winemaking are included in the terroir concept, much of the past research into terroir and terroir zoning (e.g., 1219) has focussed solely on biophysical factors, especially climate and soil physical properties. It has also relied on analysis of entire winegrowing regions rather than focussing just on the land used for winegrape production. In many instances, especially in ‘Old World’ countries, the notion of terroir aligns closely to national appellation systems, such as the French ‘Appellation d’Origine Controllée’ (AOC) system (see Robinson and Harding [20] for a summary) and to the specification of product characteristics—as in the ‘Protected Designation of Origin’ (PDO) system used more broadly throughout the European Union [21]. As a consequence, terroir zoning research in these regions has generally been constrained by conformity to the boundaries of existing denominations. This, coupled with a reliance on heuristics/presumptions of inter-regional difference and expert opinion of wines, in addition to land classification approaches based on thematic mapping [22], has led to the distinctiveness of some terroir zones being called into question [23, 24]. It has also led to a call for “unbiased scientific approaches” to be brought to bear on the study of terroir [25], focussed on process-based understanding of the complex functional relationships between terroir factors and the attributes of wine [26].

Recent research has demonstrated how a data-driven approach might be applied to the biophysical aspects of terroir zoning at the regional scale. This work has used methods of spatial analysis which have become common in the study of vineyard variability at the within-vineyard scale (e.g., 27, 28) and which have underpinned the development of precision viticulture [29]. Thus, in studies of biophysical variation in the Margaret River [30] and Barossa Zone [24] geographical indications (GIs) of Australia, k-means clustering of map layers, describing regional-scale variation in viticulturally important climate indices and soil properties, enabled the delineation of ‘zones’ within these winegrowing regions. These zones were suggested as an appropriate underpinning basis for subsequent sensory and chemical analysis of the wines produced in them, potentially leading to the delineation of subregions within these GIs for which ‘distinctiveness’ might be demonstrated. There is much interest in such subregionalisation in Australia, in the belief that it might convey marketing advantages to wine producers in these regions. However, it has also been suggested [24, 30, 31] that greater benefit might accrue from better understanding the various biophysical factors which affect final wines in terms of opportunities for improved management of the grape and wine production process. A key aspect of this recent Australian work was the observation that a different zonation resulted when soil and climate data pertaining to just that land which is used for winegrape production was included in the analysis, compared to when it was undertaken for the winegrowing region as a whole (i.e., the entire GI). This was important given that, in the Margaret River and Barossa GIs, only approximately 3 and 11% of the land is under vine. Note however, that neither of these studies included vineyard performance metrics or sensory or chemical analysis of wines.

In a third study conducted in the Marlborough region of New Zealand, Bramley et al. [22] collected data on the yield and harvest date of Sauvignon Blanc from approximately 525–750 vineyards over five vintages (2014–18) and used these to interpolate regional-scale maps of yield and harvest date variation. A key motivation for this work was to see whether vineyard performance in one location might be used to inform decisions in another. The seasonal maps showed remarkably similar patterns of variation in both yield and harvest date, despite interannual variation in the mean yield resulting from seasonal variation in climatic conditions. This similarity in patterns of variation was strongly suggestive of it reflecting a regional terroir, with both soil and temperature variation nominated as possible drivers of the variation in vineyard performance. However, this Marlborough study did not draw on any data describing soils or climate variation as was done in the Australian studies. Thus, the objectives of the present study were to incorporate such soil and climate data, along with the vineyard performance metrics, into a more holistic analysis of the Marlborough terroir. In particular, we wished to see whether the variation in vineyard performance could be explained by variation in biophysical factors. We also wished to take the opportunity to enhance the robustness of the previous analysis [22] through the incorporation of data from an additional vintage season and, across all seasons studied, from additional vineyards.

2. Methods and Materials

Marlborough, one of 16 local government regions in New Zealand, is located in the northeast of the South Island, with the District Council based in Blenheim, the largest town in the region (Figure 1). Whilst Marlborough is well known for its wine production, the rapid increase in the vineyard area over the last 30 years means there is currently no separately defined ‘wine region’ as in the case, for example, of the Australian GI or French AOC systems. Accordingly, for the purposes of this study, a map coverage of land under vineyard was obtained from the Marlborough District Council and from this, a regional grapegrowing ‘boundary’ was defined (Figure 1). In turn, this boundary was used as the basis for developing a 1 ha raster grid (i.e., pixels of 100 m × 100 m) which was used as the base for all subsequent mapping.

2.1. Vineyard Performance

The methods used for analysis of vineyard performance in the present study are exactly as described by Bramley et al. [22] with the exception that, here we include data for vintage 2019 in addition to 2014–18, and for all years, include data for vineyards additional to those canvassed in the earlier work. Various indices of vineyard performance were acquired, by request and in confidence, from local grape growers and wine companies, with the data collected being those which these entities routinely collect for the purposes of yield estimation, harvest record keeping, and payments to growers. Here, we again focus on yield and harvest date, the attributes for which the greatest amount of data was available. For mapping, reported harvest dates (Hrep) were converted to Julian day numbers (where 1 and 32 = 1st January and 1st February). Vines in Marlborough are generally planted in rows with a North:South orientation, with 2.4 to 3.0 m between the rows and 1.8 m within the row and are trimmed to a consistent height and width of about 1.8 × 0.4 m with a lower fruiting wire at 0.9 m from the soil level. Any effects of different row spacings on the yield per unit area were removed by expressing the data as kg/m, and the effects of seasonal variation were removed by normalising all data on an annual basis to a mean of zero (μ = 0) and standard deviation of one (σ = 1); the latter normalisation was also considered useful in protecting the privacy of growers. In any season, only data from blocks planted to Sauvignon Blanc at least three years prior were included, with all data georeferenced to the centroid of the vineyard block from which they derived; that is, the coordinates of the centre of each block were used to define the location of the block from whence the data derived. Regional scale yield and harvest date maps were then interpolated onto the base 1 ha grid using local point kriging in VESPER [32] with an exponential variogram model and a data cloud of 100 data points. Over the 2014–2019 study period, the number of data points available for map interpolation ranged from 618–1083 for yield and 524–851 in the case of harvest date. Figure 1(a) shows the geographical distribution of these data for vintage 2019; the distribution of data in prior vintages was very similar, albeit with generally fewer data points in the earlier vintages.

2.2. Climate and Grapevine Phenology

New Zealand does not have a freely available national gridded climate database such as the one used in the previous Australian work (e.g., [24]). Accordingly, the weather research and forecasting (WRF) model described by Skamarock et al. [33], as used previously in Marlborough by Sturman et al. [34], was run to simulate weather on a daily timestep at a resolution of 1 km2 for the 2013–2019 period; the additional year at the start of the study period was included to ensure the full season leading to vintage 2014. The WRF model has been successfully used and validated in many previous regional climate studies, including other vineyard regions [35]. Validation of WRF model simulations relevant to viticulture has also been undertaken specifically for the Marlborough region [3638], and the results used to correct any bias in model predictions. The variables modelled here were daily rainfall, from which growing season rainfall (GSR) was calculated, and daily temperature, from which the mean growing season temperature (GST) and season-growing degree days (GDD; base of 10°C) were calculated, with ‘season’ notionally defined as September to April. Note however, that to facilitate alignment to phenological modelling (see below) we used a season start date of 29th August rather than 1st September; 30th April was used as the season end date.

The daily temperature data generated above were used as input to phenological modelling. The dates of flowering (DOF) and of veraison (DOVN) were modelled using the grapevine flowering veraison model of Parker et al. [39] with values of for flowering and veraison of 1282 and 2528, respectively [40]; this model has been shown to perform well in characterising the phenology of Marlborough Sauvignon Blanc [4143]. The estimated date of harvest (Hest), defined for this purpose as the date on which the fruit reached a total soluble solids (TSS) content of 200 g/L, was calculated using the grapevine sugar ripeness model of Parker et al. [44]. A 5 × 5 pixel bilinear smoothing was applied to the 1 km2 outputs from both the WRF and phenological models, and these were then resampled to the 1 ha base grid for mapping alongside the other map layers developed. In addition, using simple map algebra in ArcGIS (v. 10.7.1; ESRI, Redlands, CA, USA), we also calculated the duration of the period between flowering and veraison, denoted here as ‘Growth,’ between veraison and harvest (Hest), denoted here as ‘Ripen,’ and a ‘harvest error’ (HarvErr), the difference between Hrep and Hest.

2.3. Soils

Soil property data were provided by New Zealand Landcare Research—Manaaki Whenua. Despite some recent activity in digital soil mapping in New Zealand [45], the available data for the Marlborough region (https://smap.landcareresearch.co.nz/) derived from conventional reconnaissance soil survey conducted in the 1980s and 90s. Much of this mapping makes use of categorical, rather than numerical data, and is only available in polygon (i.e., shapefile) rather than raster format. Nonetheless, useful and useable soil data were accessible for soil texture—the contents of sand, silt, and clay in the <2 mm fraction—and also for stone content (i.e., >2 mm). Because soil survey is uncertain, inasmuch that soils are variable over short distances, soil survey in New Zealand has employed a “soil family” and “soil sibling” approach [46] such that each mapped polygon (i.e., each soil mapping unit) may comprise more than one sibling; as such, the published maps have a probabilistic element to them as is also common, for example, in Australian reconnaissance soil survey [31, 47]. For the present study, we assumed that the soil properties of the dominant soil sibling in each polygon were those of the entire polygon.

The data for soil texture were provided for each ‘functional horizon’ with the depth of these horizons also reported. Accordingly, and notwithstanding the depth basis of soil hydrological properties (see below), we calculated profile weighted mean values for these measures of soil particle size to a maximum depth of 80 cm using the depth and soil property values for each functional horizon. In the case of soils that were deeper than 80 cm, we assumed that the functional horizon which coincided with a depth of 80 cm only reached that depth and that the mean soil property value reported for that functional horizon was appropriate to it being no deeper than 80 cm. For soils shallower than 80 cm, the profile weighted mean soil property values were calculated to the maximum depth of the deepest functional horizon. Of note in this regard is that “mean rooting depth” was also reported, with a sizeable proportion of the soils being listed as having rooting depths of >100 cm. However, most of the soil survey work conducted in Marlborough was done before winegrape production became a dominant land-use in the region. This is important given that the development of Marlborough as a winegrowing region involved minimal land re-forming, and that rooting depth is likely to be crop-specific and is also likely to be affected by the use of irrigation. Anecdotal evidence supports the view that in most locations in Marlborough, irrespective of soil depth, the majority of grapevine roots occur within the top 80 cm of the soil profile [48]. Thus, we calculated profile weighted means over this depth range.

In addition to soil texture, data were also available for the available water capacity (AWC—the amount of water potentially available for plant growth that can be stored in the soil) to 30 and 60 cm depth, with an estimate of profile available water (PAW—i.e., AWC to 100 cm depth or to a physical root barrier if one was present at less than 100 cm) also available. In each case, these were defined as the difference between water holding capacity (%) at −10 kPa and −1500 kPa in the included functional horizons, weighted by their thicknesses. As such, they are consistent with our profile weighted value for soil texture, notwithstanding the different depth ranges. The estimated AWC data used here derived from the approach of McNeil et al. [49] and accounted for the presence of stones. In addition, the categorical “soil drainage status” was available, classified in terms of very poorly, poorly, imperfectly, moderately well, or well drained [50].

Recent soil survey activity in Marlborough (G. Grealish—pers. comm.) suggests that the line work (i.e., polygon boundaries) in the existing soil mapping remains accurate. We therefore assigned the polygon-based soil values for each soil property to each coincident pixel in the base raster which aligned with that polygon. However, as can be seen in Figure 1(b), the area for which soil data were available was somewhat smaller than the area under vine. Accordingly, our soil-based analysis was confined to those parts of the district for which coincident soil and vineyard data were available.

With the exception of drainage status, the soil data were provided as numerical values. However, it was apparent from their distributions that the data were not continuous, a problem that was likely compounded by our method of assigning individual values derived from a polygon of many ha in area to all the coincident 1 ha pixels which aligned with it. In other words, large numbers of pixels could have the same soil property values. Furthermore, for some soil properties (especially the content of stones and silt), there were many pixels containing values of “zero.” For these reasons, for the purposes of clustering the various map layers (see below), the soil data for all soil properties were converted to normal scores [51] prior to cluster analysis. This uses a ranking process to calculate standard normal quantiles of the same size as the original data set. To avoid the introduction of bias in this process, prior to ranking each soil property, the order of all pixels in the dataset was randomised. The results of the cluster analysis were then interpreted in terms of real values by using the normal scores data as a ‘lookup’ table.

2.4. Topography

Elevation data were acquired at 1 m resolution by Land Information New Zealand—Toitu Te Whenua (LINZ) using airborne LiDAR in 2014, 2018, and 2020 (https://data.linz.govt.nz/search/?q=Marlborough+lidar) with each dataset covering different parts of the region, albeit with some overlap. These data were used to create a single digital elevation model (DEM) by mosaicing them in ArcGIS (v. 10.7.1; ESRI, Redlands, CA, USA). However, the LINZ dataset did not cover the upstream parts of the grapegrowing area in the Awatere Valley, nor much of the hills separating the Wairau and Awatere valleys. Accordingly, these areas were in-filled using the 8 m resolution DEM available from the Marlborough Regional Council which had been interpolated from a dataset of 20 m contours (https://data.linz.govt.nz/layer/51768-nz-8m-digital-elevation-model-2012/).

2.5. Spatial Analysis

From the above, we had a dataset of vineyard performance data (yield, Hrep), along with data for climate indices (GSR, GST, and GDD) and estimates of the date of key phenological stages (DOF, DOVN, and Hest) and derived phenological data (Growth, Ripen, and HarvErr) across the 2014–19 vintage period. Also available were data for soil attributes (texture, drainage status, AWC, and PAW) and elevation. These data were all either interpolated or sampled to the same 1 ha base raster grid except in the case of soil data for which the areal extent of the raster was matched to the area of data availability. Nonetheless, all map layers had identical alignment. This enabled similarity in patterns of spatial variation amongst these properties and across vintages to be examined using k-means clustering in JMP (v.16.0.0, SAS Institute Inc., Cary, NC, USA) with the number of clusters allowed to vary from two to five. The optimum number of clusters was selected using the cubic clustering criterion [52] and when the optimum number was initially identified as five, the analysis was rerun to larger cluster numbers to enable an unconstrained optimum to be identified. All other spatial analyses, along with map display, were done using the ArcGIS software suite (v. 10.7.1; ESRI, Redlands, CA, USA).

3. Results

3.1. Variation in Vineyard Performance

Consistent with the previous results [22], patterns of variation in both yield (Figure 2) and harvest date (Figure 3) were remarkably stable over the six years of the study, despite the observed differences between their viticultural seasons (Table 1). Thus, some parts of Marlborough can be seen to be inherently lower or higher yielding than others with a notable difference, on average, between the Wairau and Awatere valleys. In both valleys, upstream areas tend to be lower yielding than more central and downstream areas, and a marked band of higher yielding vineyards running NW-SE across the lower Wairau is consistently seen (Figure 2). Similarly, harvest dates in the Awatere Valley are generally later than in the Wairau Valley, with harvest in the Central Wairau Valley, the oldest established winegrowing area, occurring the earliest (Figure 3).

Clustering the map layers for yield and harvest date (Figure 4) emphasises the noted within-region variation and further suggests a distinction between the two valleys. In the case of yield, the optimal number of clusters based on the CCC was four (Figure 4(b)) even though the three-cluster solution, which is very similar to that reported by Bramley et al. [22], shows a more consistent rank order of the clusters (Figure 4(a)). Similarly, in the case of the clustering of harvest dates, where the two-cluster solution (Figure 4(c)) was identified as optimal based on the CCC, the three-cluster solution (Figure 4(d)) is meaningful inasmuch that the rank order of the cluster means is consistent across the six years of the study. These results, along with those for when the yield and harvest date maps are clustered together (Figure 4(e)), support the view that, in general, the Central Wairau Valley, which comprises both lower and higher yielding vineyards, is always harvested earliest, whilst most of the Awatere Valley, along with the upstream parts of the Wairau and tributary valleys are harvested last—even though these are also the lowest yielding. This observation suggests the possible importance of temperature to these patterns of regional variation. On the other hand, if temperatures (and incident sunlight) were constant, one might expect higher yielding areas to ripen later and so be harvested later than lower yielding areas. Figure 4 does not otherwise suggest a clear interaction between yield and harvest date (discussed further below).

Despite the greatly enhanced dataset which underpins the various maps in Figures 24 compared to the previous analysis [22], the interaction between the density of data and use of local kriging means that the map confidence intervals do not allow us to identify the statistical significance of the difference between the cluster means based on the median kriging variance [53]. Accordingly, to examine the significance of between-cluster differences, the raw point data for yield and harvest date (e.g. Figure 1(a)) were overlaid on the results of the cluster analysis (Figure 4), and the raw data were reanalysed to test for the statistical significance of differences in the means (Tukey–Kramer test) of the identified clusters. The results are presented in Tables 24, with a colour coding to assist with matching to the map legends. Note that for this analysis, we only used data from vineyards for which we had both yield and harvest date information; we also used the actual yields (kg/m) rather than normalised data.

Comparison of Table 2 and Figures 4(a) and 4(b) suggests that the different yield clusters identified do reflect the generally statistically significant yield differences in the raw data derived from vineyards located in each cluster. A very similar conclusion may be drawn in terms of the harvest date (Table 3 and Figures 4(b) and 4(c)). When a similar analysis is done for the combined yield and harvest date clusters, the results (Table 4) again mirror those from the cluster analysis (Figure 4(e)) and, as such, reflect the lack of interaction between the yield and harvest date over the region as a whole; that is, different factors appear to drive yield variation compared to harvest date variation. Of note is the fact that, when a similar analysis is done to test for differences between the Wairau and Awatere valleys (i.e., ignoring the cluster analysis), the two valleys present as clearly different (Table 5).

3.2. Variation in Climate

Consistent with temporal stability in patterns of variation in vineyard performance (Figures 24), patterns of spatial variation in both temperature (Figure 5) and rainfall (Figure 6) were also consistent across the six years of the study, even though the absolute values varied from year to year. This was the expected result given the strong reliance of the WRF model at the local level on topographic variation and distance from the sea, coupled to the effects of prevailing weather systems. This is also the reason why we only show GDD in Figure 5 since its patterns of spatial variation were essentially identical to those for GST. Of note is the apparent difference between the patterns of variation in annual rainfall (Figure 6) and those for vineyard performance (Figures 24), in spite of the general observation within both the Wairau and Awatere valleys of a strong S-N rainfall gradient (Figure 6; in the Wairau Valley, rainfall roughly doubles between Blenheim and the north bank of the Wairau River, 10 km due north), and the increase in rainfall with distance upstream (approximately 16 mm/km in the Wairau Valley). On average, the Awatere Valley does not present as markedly cooler than the Wairau Valley (Figure 5) suggesting that, contrary to the comments above, some factors other than temperature accumulation might be responsible for the lower yields in the Awatere Valley.

3.3. Variation in Phenology

Patterns of variation in grapevine phenology closely followed those of temperature. This was the expected result given the dependence of the phenological models of Parker et al. [39, 40, 44] on temperature. Thus, Figure 7 shows mean predicted dates of flowering, veraison, and harvest (i.e., TSS of 200 g/L; Hest) across the six study years, whilst Table 6 provides an indication of the interannual variation. Of note is that, whilst the patterns of spatial variation are distinct, the range of variation in DOF, DOVN, and Hest is narrow in any given year (Table 6). Furthermore, whereas the pattern of spatial variation in the duration of the flowering to veraison growth period (Growth; Figure 7(d)) aligns closely to that of heat accumulation during the whole growing season (GDD; Figure 5), the pattern of spatial variation in the duration of the shorter ‘Ripen’ period between veraison and Hest, the date at which fruit reach a TSS of 200 g/L (Figure 7(e)), shows a less organised pattern of variation, in spite of the general temperature dependence of vine phenology (Figures 7(a)–7(c); [39, 44]). This is likely a consequence of the difference in length between the short ‘Ripen’ and much longer ‘Growth’ phenophases, their relative proportions of the length of the overall growing season, and heat summation over the whole season. However, of greater interest in the context of terroir is that the largest values of HarvErr (i.e., Hrep-Hest; Figure 7(f)) occur in the Awatere Valley, and in the high yielding strip in the lower Wairau (Figures 2 and 4); that is, in the Wairau Valley, and contrary to the suggestion (above) of no interaction between the yield and harvest date, the delay in harvest beyond a TSS of 200 g/L appears related to higher yields. Since the patterns of variation in Hest across the two valleys (Figure 7(c)) are similar to those of temperature (Figure 5), the delayed harvest in the lower yielding Awatere Valley (Figure 7(f)) cannot be attributed to temperature, as was suggested by Bramley et al. [22].

3.4. Soils

Patterns of variation in available water capacity in Marlborough soils are very similar whether expressed to a depth of 30, 60, or 100 cm (Figures 8(a)–8(c)) and show an inverse relationship to sand content (Figure 8(e)) and stoniness (Figure 8(h)); that is, as is expected, soils with higher contents of sand and stones, such as predominate upstream of the central Wairau Valley, have lower AWC. These soils are also well drained (Figure 8(d)). Conversely, and as is expected in a floodplain, lower Wairau soils have higher contents of silt (Figure 8(f)) and clay (Figure 8(g)) and these areas are also less well drained (Figure 8(d)). An exception occurs in the northeasternmost area of Marlborough around Rarangi where the soils reflect relic beaches and are characterised locally as being composed of ‘pea gravel.’

Consistent with the above, clustering of soil texture data (Figure 9(a)) divided the region into two clusters of soils with either few stones, comparatively low sand contents, and higher silt and clay contents on the one hand, and much sandier, stonier soils on the other. When the soil available water data were clustered (Figure 9(b)), four clusters were identified with a consistent rank order of AWC amongst the clusters irrespective of the soil profile depth increment. When all the soil properties were clustered (not shown), an almost identical delineation resulted as the more parsimonious analysis shown in Figure 9(c) which included just PAW (i.e., AWC to 1 m depth) and the profile weighted contents of silt and stones to a depth of 80 cm. Comparison of this result with Figure 8(c) suggests that the variation in the available soil properties can be characterised by variation in PAW—which is also the expected result given the dependence of PAW on texture.

3.5. Integrating Variation in Vineyard Performance with Soil and Climate Variation

Clustering a map of mean yield (i.e., the mean of the six maps shown in Figure 2) with PAW using k-means splits the region into two, based on CCC; an above-average yielding area with a mean PAW of 176.6 mm and an area of approximately average district yield (μ ≈ 0) with PAW of 98.3 mm (Figure 10(a)). Repeating this analysis with the yield maps for individual years (Figure 2) rather than the mean, gave a very similar result, as did using silt content as the soil property included in the analysis, instead of PAW (not shown). However, when the yield maps for individual years were clustered with PAW and the contents of both silt and stones, no maximum cluster number was reached; that is, even when the clustering was allowed to run to 20 clusters, 20 was the identified optimum number based on the CCC. A possible explanation may be the general trend for yield to increase with increasing PAW (Figure 10(b)). Nonetheless, the results are consistent with the apparent alignment of the pattern of variation in both drainage status (Figure 8(d)) and silt content (Figure 8(f)) with that for yield (Figures 2 and 4). Thus, the higher yielding area in the lower Wairau aligns closely with these poorly drained soils with high silt content—an instance where, contrary to what in other regions might be ‘conventional wisdom,’ poorer drainage promotes higher vine vigour and yield.

Clustering GSR with either yield, harvest date, or both yield and harvest date together did not suggest any impact of regional rainfall distribution on vineyard performance. In both the Wairau and Awatere Valleys, the wettest areas in the upstream parts of the catchments tend to be the lowest yielding and are also harvested later in the season. Including GSR in the cluster analyses otherwise simply tended to reflect the patterns seen in Figure 6 without clear interaction with vineyard performance. Thus, the lower yields and later harvests seen in upstream areas are unlikely to be caused by their wetness and are likely driven more by their lower temperatures (Figure 5). Likewise, including the duration of the modelled ripening period (Ripen) had no effect in separating clusters. In contrast, including HarvErr (the difference between the modelled and actual harvest date) in the cluster analysis did align meaningfully with the patterns seen in many of the other maps (Figure 11). Thus, clustering the mean yield and HarvErr over the 6 seasons, either with or without PAW (Figures 11(a) and 11(b)) delineated two clusters with a familiar pattern. The higher yielding area in the lower Wairau is seen to also have higher values of HarvErr, but because parts of the lower yielding Awatere Valley also align with this higher HarvErr cluster, the mean cluster yield is seen to be lower than in Figures 4(a) and 4(b) or 10(a). When PAW is included in the analysis (Figure 11(b)), the mean yield of the higher yielding, greater HarvErr, and PAW cluster increases by comparison with Figure 11(a), as the Awatere area is no longer included in the analysis given the lack of soil information in that area. When the analysis is expanded to the individual years, three clusters are identified which more explicitly separate the Awatere from the Wairau, in particular highlighting the delay to harvest in the Awatere Valley (Figure 11(c)) and the association between higher yield and PAW and the likely impact of this higher yield in terms of later harvest. Of note is that PAW is not a discriminator between the clusters other than in regard to the high yielding strip in the lower Wairau (Figure 11(d)). Expanding this analysis to also include temperature (GDD) and harvest date (Hrep; Figure 12) provides very similar results, albeit perhaps suggesting that silt content may be a better discriminator between the soils of the different clusters (Figure 12(b), Table 7). Since the patterns seen in Figures 12(b) and 12(c) are essentially the same, it seems clear that soil variation within the Marlborough region is not a primary driver of variation in vineyard performance, other than in terms of the impact of the high PAW and silt / low sand and stone soils of the lower Wairau. Excluding the Awatere Valley from the analysis (Figures 12(d)–12(f)) lends weight to this conclusion, with comparison between Figures 12(b) and 12(e) and between 12(c) and 12(f) highlighting the distinction between the Wairau and Awatere Valleys.

4. Discussion

Much previous terroir research has focused on the delineation of so-called “homogeneous” terroir zones (e.g., [13, 18, 26]). However, just as research into vineyard variability and precision viticulture [54] has suggested that there is no such thing as a uniform vineyard, comparison of Figures 2, 3, 5, 6, and 8 with Figures 4 and 1012 supports the view that, similarly, homogenous terroir units do not exist. Recognising the pragmatic need, across scales, to organise and classify variable data in a manner consistent with its intended use, we nevertheless consider the description of zones at any scale as “homogenous” (i.e., invariant) as unhelpful, especially if a part of the objective of the classification is to support process-based understanding of the complex functional relationships between terroir factors and the attributes of wine [26]; such relationships will, of course, be subject to error. Using a technique such as k-means clustering, it is certainly possible to identify clusters or subregions in which the range of variation within the clusters is substantially less than the variation in the region as a whole. This is especially the case in a region like Marlborough, given the observed differences between the Wairau and Awatere valleys (Figures 4, 11, and 12). Similarly, Bramley and Ouzman [24] demonstrated how the soils and climate of the Barossa and Eden Valleys were different, in spite of being part of the same Barossa Zone GI. The important difference between these New Zealand and Australian studies is in the incorporation of vineyard performance metrics in the analysis—something that was not possible in the Barossa example. Of course, neither study was able to incorporate chemical or sensory analysis of wines and it is to be hoped that future work will enable this so that the implications of biophysical variation for final wines might be better understood and relevant functional relationships [26] developed. Nonetheless, both studies speak to variation in terroir at the subregional scale, with the inclusion of vineyard performance metrics lending weight to consideration of the importance of observed biophysical variation in the landscapes in which grapes are grown and wine is made.

Aside from short-range variation, an obvious reason for heterogeneity within the identified terroir zones is variation in the specific production objectives associated with different vineyard blocks and resultant variation in grower management practices, such as trellis design, pruning to particular bud numbers, crop thinning, management of disease risk, irrigation, and the timing of decisions associated with these things. Timing of harvest and its interaction with winery logistics is also a potentially large source of confounding error, as might be vine age. One of the original motivations of this study [22] was to understand whether estimates of yield made in one location could be used to inform estimates needed in other locations, given the practicalities of deploying sensors that might assist with such estimation and/or associated labour to many different vineyards at optimal times. In spite of the constraints imposed by the variation in the nature and timing of grower and winemaker management and the other factors noted above, the fact that a marked and consistent spatial structure in the regional-scale variation in vineyard performance can be noted and interpreted in the context of variation in the winegrowing environment is important. It also lends weight to the idea that, given similarities in vineyard characteristics and management, estimates of yield in one vineyard could be useful estimators of yield in other vineyards in the same zone, just as at the within-vineyard scale, zone-based sampling may be useful [29, 55]. Understanding of such variation could also be used to inform winery logistics and associated harvesting decisions, especially when coupled to understanding of variation in fruit quality [56].

Variation in harvest date (Figure 3) is largely a result of variation in temperature-driven vine phenology (Figures 5 and 7). Thus, GDD decreases as one moves upstream in both the Wairau and Awatere valleys and into the southern valleys of the Wairau Plain. Overall, the Awatere Valley is cooler than the Wairau, inasmuch that the warmer areas comprise a proportionally smaller faction of the grapegrowing area in the Awatere than is the case in the Wairau (Figure 5). However, the differences are arguably not large enough to explain the gross differences in vineyard performance between the two valleys. It is possible that some of the difference can be explained by differences in daily wind run and maximum wind speed observed between the valleys. While the Marlborough region is protected from southerly winds by the inland and seaward Kaikoura Mountains (maximum altitude 2885 m), the western Awatere valley has greater exposure. The Wairau Valley is protected by the Black Birch range of hills (maximum altitude 1500 m). The result is the average daily wind run and maximum daily wind speed are consistently 66 km/day and 5 km/h greater in the Awatere than the Wairau valleys (M. Trought–pers. comm.; https://www.mrc.org.nz/blenheim-weather-station, https://www.mrc.org.nz/awatere-weather-station), while the predominant wind direction in the Wairau Valley is westerly (at right angles to the canopies), and in the Awatere Valley it is mainly north-westerly or south-easterly, and therefore parallel to the row orientation [57]. Wind flow over grapevine canopies has received little previous attention. While wind direction across rows results in flow similar to a uniform canopy, spatial variability increases as the direction changes to become parallel to the rows, which generates greater turbulent kinetic energy and turbulent flux [58]. Similarly, the effect of wind on grapevines is poorly understood. Using artificial wind breaks, Dry et al. [59] and Bettiga et al. [60] reported increases in Cabernet Franc and Chardonnay yields, largely as a result of better budbreak and heavier bunches under sheltered conditions. Perhaps more wind in the Awatere is a part of the explanation for both its lower yield and later harvest dates? Conversely, variation in yield is associated with soil characteristics, at least in respect of the high yielding band of poorly drained, silty soils in the lower Wairau with higher PAW than in the rest of the region. Trought et al. [61] and Bramley et al. [62] have noted the impact of soil variation at the within-vineyard scale, with narrow silty hollows—a relic of the active Wairau floodplain—promoting higher vigour, but not higher yield [63], with this soil variation also having profound implications for grape and wine quality [56, 61]. At regional scale, the impact of texture was much less clear [64], with neither yield nor location within the Wairau Valley impacting on wine sensory properties in spite of differences in soil texture; fruit from an Awatere Valley vineyard had a higher methoxypyrazine concentration and herbaceous characteristics than the Wairau wines, when fruit was harvested at the same soluble solids. Conversely, Jouanneau et al. [65] analysed wine aroma compounds in “research-scale” wines made from juices of variable soluble solids collected from seven predetermined subregions within Marlborough and noted a lower methoxypyrazine concentration in Awatere wines. Whilst they noted some subregional differences in wine chemistry, no attempt was made to relate these to fruit ripeness, which influences both thiol and methoxypyrazine concentrations [66], soil properties, or other biophysical attributes. However, the basis for the subregional delineation used by Jouanneau et al. [65] is not clear and is not supported by either the present analysis or the results of Trought et al. [64], which may explain why the distinctiveness of the subregional wines was equivocal. Similarly, the anecdotal local suggestion that Marlborough be divided into the Southern Valleys, Awatere and Wairau Valleys, which seemingly derives largely from the history of Marlborough’s development, does not appear to be otherwise underpinned by data, aside from the presence of more clayey soils in the Southern Valleys (Figure 8(g)). It is also a fact that the original plantings in the Southern Valleys used an E-W row orientation in contrast to the more common N-S orientation in the remainder of the region. However, while the impact of soil properties on higher yields in the lower Wairau is clear, the previous studies, along with the available soil data (Figure 8), have generally supported the view that soil properties are not a major driver of regional-scale terroir variation in Marlborough. Arguably, this might be due to the geologically young and relatively undifferentiated soils in the region—clay contents are generally low throughout (Figure 8(g))—coupled with the need for irrigation to support commercial viticulture.

Comparison of predicted (Hest) and actual (Hrep) harvest dates (i.e., HarvErr) adds some interesting observations. Hest (Figure 7(c)) is temperature-driven and reflects the date of flowering (Figure 7(a)). In contrast, Hrep (Figure 3) is more variable with the higher yielding strip in the lower Wairau being harvested relatively late, as indicated by values of HarvErr. Of course, actual harvest dates are influenced by factors other than temperature, such as vineyard management practices. Thus, the time from flowering to veraison is increased as yields increase [67] and the time from veraison to target soluble solids is strongly influenced by yield [68], although presumably due to variation in management practices, this is not evident in comparison of Figure 7(e) with Figures 2 and 4(b). Furthermore, notwithstanding winery requirements for fruit that is fit for intended end-use [66], many harvest management decisions will be influenced by the proximity of the harvester, particularly in wetter vintages when disease risk is an important determinant of quality (A Naylor, Pernod Ricard NZ—pers. comm.). Proximity of the harvester could arguably be one reason for the generally later harvests in the Awatere Valley compared to the Wairau Valley where the majority of the wineries are located. Finally, it is interesting to note that yields in the Rapaura area of the Central Wairau are lower than adjacent areas with similar soils. This area was the first to be planted within the Marlborough region and so contains the oldest vines with potentially greater numbers succumbing to trunk disease [69] than in other parts of the region.

Previous work has led to differing conclusions as to the importance of elevation and topographic variation to terroir at different scales. In the Margaret River region of Australia, elevation did not contribute usefully to delineation of subregions within the GI [30] because Margaret River vineyards neither occur at locally characteristic elevations nor have locally characteristic slopes or aspects. Similarly, Biss [70] argued that the impact of topography on wine quality in the Chablis region of France was equivocal, yet in more mountainous areas such as the Italian Tyrol, topography may clearly have a major impact [71]. Likewise, in the much less mountainous Barossa Zone GI, topographic variation was seen to play an important role in both delineating the Barossa and Eden Valleys, and also in explaining some of the subregional differentiation within the Barossa Valley [24] and, at property scale, within the Eden Valley [31]. The aforementioned impact of soil property variation at the within-vineyard scale in Marlborough on vine vigour and fruit quality is directly attributable to topographic variation as demonstrated by Bramley et al. [62] and Trought and Bramley [56]. In the Wairau Valley, the active floodplain is characterised by a pattern of silty hollows (Wairau series) dissecting the sandier, more gravelly soils (Rapaura series) that predominate [72]. These soil series represent approximately 3130 ha of the Wairau Plains, but the effects of the silty hollows are not evident at the scale at which the soil property data are available (Figure 8) and the underpinning soil survey was conducted; the 1 ha base raster used for the present mapping also presents a difficultly in this regard. The hollows run predominantly in an approximately east-west direction, while row orientation is generally north-south. Thus, the full range of within-vineyard variation is commonly expressed in a single row. These scale effects (regional vs. within-vineyard) are likely the reason for the fact that, in the present study, soil properties did not clearly impact on regional-scale variation in vineyard performance beyond the high yielding part of the lower Wairau. Much the same conclusion can be drawn in respect of topographic variation; there is a mismatch between the high resolution (1 m) LINZ DEM and the resolution of other available data which is why, aside from the hills which enclose the Wairau and Awatere Valleys, topography is not evidently a strong driver of regional terroir variation (Figure 13). Further research aimed at understanding how to integrate terroir expression at different scales would therefore be valuable, especially in respect of topography and soils.

Finally, a key reason for Australian interest in subregionalisation is the belief that it may promote marketing advantages to wine producers in different parts of Australia’s (generally large) GIs through the ability to demonstrate the ‘distinctiveness’ of their wines and so use their terroir as the basis of the ‘story’ used to sell them. However, Charters et al. [73] cautioned against the difficulty of promoting “territorial” brands when, understandably, most producers attach primacy to their proprietary brands. Organisations with oversight of the territorial brand are also generally distinct from individual producers. The observation of differences between the Wairau and Awatere Valleys is therefore interesting given that ‘Marlborough Sauvignon Blanc’ has developed a significant international presence in the wine trade. In other words, for reasons of marketing, Marlborough producers might not wish to pursue subregionalisation, even though Table 5 and Figures 4, 11(c), 12(b) and 12(c) present a strong justification as to how they could. Conversely, Pinu et al. [74] used a metabolomics approach to assess and compare 400 Sauvignon Blanc juice samples from around New Zealand, 75% of which came from Marlborough. They concluded that seasonal variation was of greater importance than geography in discriminating the characteristics of the samples, although they did not present a within-Marlborough analysis of the Marlborough samples. Nonetheless, perhaps Marlborough presents a case where, contrary to Charters et al. [73], the territorial brand is stronger than, or at least as strong as, its proprietary brands with consistency of style achieved by blending wines from various subregions of Marlborough. If so, the present study lends weight to the idea [25, 31] that understanding terroir has more to offer the optimisation of grape and wine production systems than to being used for marketing objectives.

5. Conclusions

Analytical techniques that have previously been applied to studies of within-vineyard variability and the development of precision viticulture are valuable tools in assessing regional-scale variation in biophysical variation and vineyard performance, which might impact on and reflect a regional/subregional terroir. Their use in the present study strongly suggests that the Marlborough region has a characteristically variable Sauvignon Blanc production with crop phenology and harvest date strongly influenced by variation in temperature, and yield variation impacted by soil properties, albeit less distinctly than is apparent at the within-vineyard scale. A key part of this is the apparent distinction between the Wairau and Awatere Valleys which, hitherto, have been considered together as parts of a single Marlborough region. The results from this study, which to the knowledge of the authors is the first quantitative integration of vineyard performance and biophysical metrics as a means of evaluating terroir, has potentially important implications for the management of both vineyard operations and winery logistics, for wine marketing and potentially, for whole-of-industry planning around expansion or contraction. It also lends weight to the idea that estimates of vineyard performance in some parts of the region may be used to predict performance in others. Accordingly, a coordinated collection of vineyard performance metrics is encouraged for all regions to better understand their terroir—especially as most grapegrowing businesses collect such data as a matter of course.

Data Availability

The vineyard performance data used in this study are confidential and are not available. Details of the soil and climate data used are provided within the article.

Disclosure

The present address of G. J. Grealish is CSIRO, Black Mountain Science and Innovation Park, Canberra, ACT 2601, Australia.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was jointly funded by CSIRO, the University of Canterbury, Manaaki Whenua—Landcare Research, New Zealand Plant and Food Research, Innovative Winegrowing, and the New Zealand Ministry of Business Innovation and Employment (MBIE) through an ‘Endeavour grant’ provided to Lincoln Agritech Ltd. (MBIE contract LVLX1601). The authors are most grateful to the Marlborough District Council for providing access to their coverage of vineyard holdings, and in particular, to the following wine companies who contributed data for the work: Babich, Cloudy Bay, Constellation Brands, Delegat, Forrest, Giesen, Indevin, Pernod Ricard, Treasury Wine Estates, The Marlborough Grape Growers Cooperative, The Wine Portfolio, Villa Maria, Weta, Whitehaven, Wither Hills, and Yealands. Soil data were obtained from Manaaki Whenua—Landcare Research. Drs Amber Parker (Lincoln University, New Zealand) and Brent Sams (E&J Gallo Winery, USA) made valuable comments on a previous draft of the manuscript.