Effects of microtopography on within-field spatial variation of inherent soil properties in Andosols of Tokachi district, Hokkaido

ABSTRACT

In large-scale upland cropping systems, understanding the magnitude of soil property variations and mapping the distribution are essential to improve soil nutrient management. The Tokachi district of Hokkaido has large-scale upland fields, and the soils are dominated by Andosols that inherently have characteristic clay minerals and high soil organic carbon content. In the region, microtopography is known to exist that can affect soil formation and within-field variation of the inherent soil properties. However, limited information exists about the relationship between microtopography and the spatial variation of inherent soil properties that are important for nutrient dynamics. Therefore, this study aimed to assess the effect of microtopography on the within-field variation of inherent soil properties in Andosols of Tokachi district, Hokkaido. Within a 6.6-ha field with an elevation difference of 3.2 m, four soil profiles were surveyed at different topographical positions, and 530 surface soil samples were collected at 12-m intervals. Samples were analyzed for total carbon (TC), acid-oxalate extractable Al, Fe, and Si (Al_o, Fe_o, and Si_o), pyrophosphate extractable Al and Fe (Al_p and Fe_p), cation exchange capacity (CEC), and phosphate absorption coefficient (PAC). Geostatistical analysis was applied to understand the spatial structure of the surface soil properties. All soil profiles were dominated by volcanic parent materials, and they were classified into different types of Andosols. An exceptional variation of inherent soil properties was found in the surface soil, with TC content ranging from 35.2 to 124 g kg⁻¹ and CEC from 20.3 to 55.7 cmol_c kg⁻¹. Geostatistical analysis showed that TC content varied over a large distance (range: 226 m) while allophane content varied at a shorter distance (range: 118 m). Our findings revealed that microtopography affected the soil water dynamics and erosion of the surface soil particles which caused the variation of TC, CEC, and allophane. Further studies must be conducted to understand the effects of the variation of these inherent soil properties on soil nutrient availability and crop productivity.

KEY WORDS:

• Andosols
• inherent soil properties
• upland
• microtopography
• spatial variability

1. Introduction

Precision management of agricultural inputs is becoming more critical in arable cropping systems since the constraints over water, soil, biodiversity, and land are even more significant (Pretty and Bharucha 2014). In Japan, the increasing cost of fertilizer and the recent report on the accumulation of soil nutrients (Obara and Nakai 2004; Tani and Kinoshita 2020) are inducing farmers to optimize fertilizer application to improve profitability and reduce environmental damage. For large-scale agricultural production, understanding the source and scale of soil variability is important for the improvement of soil and crop management (Di Virgilio, Monti, and Venturi 2007).

Soil properties vary spatially from a field scale to a larger regional scale, affected by soil-forming and anthropogenic factors (Cambardella and Karlen 1999). The soil properties can be classified into inherent and dynamic, and the former ones are soil properties that do not change or change very little with management. They have a high degree of spatial variability due to the combined effect of physical, chemical, and biological processes operating in different intensities and at different scales (Chopra and Chakraborty 2008; Kilic, Kilic, and Kocyigit 2012). At the field scale, the spatial variations of the inherent soil properties are likely to affect the spatial variation of the dynamic properties including soil nutrient dynamics that are important for crop production and environmental conservation.

The Tokachi district of Hokkaido is one of the most important upland crop production areas in Japan on a large-scale farm with 45.7 ha on average, which is 21 times larger compared with the rest of the country (MInistry of Agriculture, Forestry and Fisheries, Japan 2020). The surface topography is a complex of alluvial fans and floodplains on river terraces of different formation ages (Kikuchi 2008). The distribution of parent materials depends on the formation age of the river terrace, varying from alluvial materials to volcanic ash (Hatano, Shinjo, and Yusuke 2021; Kato et al. 2019; Kikuchi 2008; Tani et al. 2010). The soil in the high terrace was mainly formed from volcanic parent materials deposited during the Pleistocene and Holocene. In the middle and low terraces, the soils are dominated by volcanic parent materials on the surface, deposited during the late Pleistocene and Holocene, overlying alluvial deposits in the subsoil. However, the thickness of the volcanic materials in the middle terrace is greater than that observed in the low terrace. In the lowland, the soils are younger than those in the terraces, dominated by alluvial materials recently deposited by existing rivers (Kikuchi 2008).

Microtopography is also known to affect the spatial variation of soil properties on a smaller scale. In the Tokachi district, many gentle slopes are found, especially in the river terraces that were affected by periglacial processes, such as ice erosion and sand dune formation (Kikuchi 2008; Niwa et al. 2011). These abundant slopes are responsible for the microtopography commonly found at the field scale. Microtopography refers to soil surface variation, including patterns of elevation at many spatial scales, encompassing both vertical and surface roughness (Larkin, Vivian-Smith, and Zedler 2006; Moser, Ahn, and Noe 2007). Imai (1964) suggested that the spatial variation of soil properties in the Tokachi district might be because the volcanic ash and alluvial deposits were eroded and redeposited during the Holocene and Late Pleistocene ages. These pedogenic processes might have caused a huge impact on the spatial distribution of inherent soil properties at a within-field scale since the field size can be as large as a few hectares or more (Niwa et al. 2011).

About 58% of the arable soils are classified as Andosols and about 28% as Lowland soils based on the Comprehensive Soil Classification System of Japan First Approximation (Kanda et al. 2016). Andosols are known to have high soil organic carbon (SOC) content that provides high water-holding capacity, soil buffering, and CEC (Matus et al. 2014). Andosols also have an abundance of poorly crystalline clay minerals such as allophane, imogolite, and ferrihydrite as well as an Al-humus complex, which determine various soil functions important for crop production along with SOC (Shoji, Dahlgren, and Nanzyo 1993). For example, they show a very high phosphate fixation capacity that limits plant-available phosphate, the major yield-limiting factor in the region (Shoji, Dahlgren, and Nanzyo 1993). In Hokkaido, there is already a fertilizer recommendation system that is soil type specific for a particular crop type, cultivar, region, and cultivation method (Department of Agriculture, Hokkaido Government 2020), which aims to determine an optimum fertilizer rate.

Within-field variation of soil properties has been assessed widely for different soil types. Early studies such as Cambardella et al. (1994) showed a spatial variation of 28 soil parameters for two fields under Mollisols with 36 and 96 ha, respectively. There were large spatial variations of inherent soil properties, including SOC and soil texture, with a coefficient of variation (CV) of about 35%, as well as for plant available nutrients (35 ≤ CV ≤ 74%). The magnitude of spatial variation was dependent on the landscape features of each field, and the within-field variation of inherent soil properties was affected by landscape positions such as summit, backslope, and depressional positions. In a more recent study in Acrisols (IUSS Working Group 2015), the spatial variation of SOC was also affected by elevation in combination with land use and parent materials, while the CEC was only affected by land use when evaluated at the regional scale (932 km²; Song et al. 2020). These types of studies are very limited for Andosols that inherently show very high SOC content and distinct mineralogical characteristics. In the Tokachi district, some attempts were made by Niwa et al. (2019), who reported soil C varying from 29.6 to 82.4 g kg⁻¹ and phosphate absorption coefficient (PAC) ranging from 1,070 to 1,750 in one field where the parent materials were volcanic ash in the topsoil and alluvial deposits in the subsoils. However, within-field variations of inherent soil properties other than soil C, such as clay mineralogy and associated soil properties important for crop production, have not been evaluated in detail for upland crop fields under Andosols. Also, the relationship between microtopography and the spatial variation of soil properties has not been assessed in the region. Therefore, this study was conducted to assess the effect of microtopography on the within-field variation of inherent soil properties in Andosols of Tokachi district, Hokkaido.

2. Materials and methods

2.1. Study site

The study was carried out in an experimental field within the Obihiro University of Agriculture and Veterinary Medicine located in Obihiro City, Tokachi district, Hokkaido prefecture, Japan (42°51’43.4“N, 143°09’57.5“E). The region is within a humid continental climate (Dfb; Peel, Finlayson, and Mcmahon 2007). The 30-year average annual precipitation and air temperature were 926 mm and 7.2°C, respectively, with a mean monthly temperature ranging from −7.1°C in January to 20.3°C in August (Japan Metrological Agency 2022). The region is surrounded by the Hidaka range in the West and the Daisetsuzan Volcanic Group in the North. Rivers have created alluvial fans and, subsequently, river terraces and floodplains in the Pleistocene. The river terraces can be classified according to the age of formation. The lowland adjacent to the river courses is the youngest and lacks volcanic materials due to recent flooding and only contains alluvial deposits, which are mostly classified into Lowland soils (Kanda et al. 2016). On the terraces, most of the soils are classified into Andosols, but the composition of the parent materials within the profiles is dependent on the formation ages of the terraces (Kikuchi 2008). The study field was located on the middle terrace (Kamisatsunai I Plain; Tokachi Research group, Tokachi Research Group, 1978), where the parent materials are known to be dominated by the late Pleistocene and Holocene volcanic materials. The field was 550 m in the N-S direction and 120 m E-W direction, a total of 6.6 ha (Figure 1a). An unmanned aerial vehicle model Phantom 4 RTK (SZ DJI Technology Co., Ltd. 2023, Shenzhen, China) equipped with a 1-inch, 20-megapixel CMOS sensor was utilized to collect the surface soil images of the field from a 25 m height, at a speed of 2.0 m/s and 80% of the photos overlap on both front and side. Afterward, the field image was used for making the digital elevation map (Figure 1b) in ArcMap 10.6.1 (Esri Inc. Berkeley, California). An undulating topography was observed within the field, with two ancient sand dunes in the northwest and southeast corners (Figure 1(b,c)), creating a 3.2-m elevation difference. The ancient sand dunes were formed with Eniwa-a (En-a) volcanic materials deposited 19,000 to 21,000 years B.P. (Machida and Arai 2003), which were secondarily sedimented by wind (Kimura et al. 1970). This is a common feature in the middle and high terraces in the central to the southern part of the Tokachi district.

Figure 1. (a) Map of surface sampling locations and soil profile survey sites, (b) elevation map, and (c) schematic diagram of soil profile survey sites in relation to microtopography. Source of the base map: Esri, Maxar, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS user community (Esri Inc. 2018).

2.2. Soil profile survey

Four soil profiles were prepared in May 2021 at different topographical positions (Figures 1(a,c)) to understand the soil formation process and especially the effect of microtopography on soil properties. Profile 2 was at the lowest elevation, while profile 4 was at the highest elevation on top of the ancient sand dune. Morphological features of the profiles were recorded in the field, and soil samples were collected from each identified horizon for mineralogical and physicochemical analysis.

2.3. Surface soil sampling

A total of 530 surface soil samples were collected in May 2020 and April 2021 (Figure 1a) using an automated soil sampler (Wintex 1000, Wintex Agro, Thisted, Denmark) mounted on an all-terrain vehicle. The autosampler was also equipped with a high-precision Global Navigation Satellite System receiver (DG-PRO1RWS; Biz Station Co., Ltd., Nagano, Japan) with a Virtual Reference Station System at 1 cm horizontal accuracy. Soil samples were collected using an equilateral triangular grid design with 12 m on each side (Figure 1a) and depths of 0–25 cm. Each soil sample was composed of 10 subsamples that were collected at every 20 cm, starting at a distance of 80 cm before the georeferenced sampling point, along a straight line in the North–South direction.

2.4. Soil analysis and data calculation

All soil samples were air-dried and passed through a 2-mm sieve, followed by the removal of visible organic residues. Total carbon content (TC) was measured using a dry combustion method with a Vario EL III analyzer (Elementar Analysensysteme, Hanau, Germany). Acid-oxalate extractable aluminum, iron, and silicon (Al_o, Fe_o, and Si_o) were extracted using a 0.2 mol L⁻¹ acid ammonium oxalate solution at pH 3.0 after shaking for 4 hours in the dark (Blakemore, Seatle, and Daly 1987). Pyrophosphate extractable aluminum and iron (Al_p and Fe_p) were extracted using a 0.1 mol L⁻¹ sodium pyrophosphate solution at pH 10.0 after shaking for 16 hours (Blakemore, Seatle, and Daly 1987). Both acid-oxalate and pyrophosphate extracted elements were measured by an inductively coupled plasma emission spectrophotometer (ICP-AES; ICPE-9820; Shimadzu Corporation, Kyoto, Japan). Allophane content was calculated according to the equation proposed by Mizota and van Reeuwijk (1989), which is widely used. The allophane content in the soil is obtained by the relation between aluminum and silicon ratio ((Al_o-A_p)/Si_o) and the silicon content in the allophane (Si_o). The values were converted to grams per kilogram for better comparability with other soil properties.

(1) $\mathrm{Allophane}(\mathrm{gkg}-1)=\frac{(100\times \mathrm{S}{\mathrm{i}}_{\mathrm{o}})}{\left[23.4-5.1\times (\mathrm{A}{\mathrm{l}}_{\mathrm{o}}-\mathrm{A}{\mathrm{l}}_{\mathrm{p}})/\mathrm{S}{\mathrm{i}}_{\mathrm{o}}\right]}$(1)

CEC was measured using the ammonium acetate extraction method at pH 7.0 (Schollenberger and Simon 1945). PAC was measured by reacting a soil sample with a phosphate retention solution (13.44 mg P₂O₅ mL⁻¹) at pH 7.0 for 24 hours, followed by determining the remaining phosphate in the solution using the molybdate yellow method (Nanzyo 1997).

2.5. Statistical and geostatistical analysis

Descriptive statistics and Pearson correlation coefficients were calculated using JMP Pro 16.1 (SAS Institute Inc., Cary, North California, SAS Intitute Inc. 2023). The gstat package (Pebesma 2004) in the R statistical computing environment (Version 1.4.1106, RStudio, PBC, RStudio Team. 2021) was used for all geostatistical modeling. Semivariance (γ(h)) is half of the variance between two locations that are separated by a distance h called lag, and m(h) is the number of pairs at the distance h, and it was calculated as equation 2:

(2) $\mathit{\gamma }\left(\mathrm{h}\right)=\frac{1}{2\mathrm{m}(\mathrm{h})}{\sum }_{\mathrm{i}=1}^{\mathrm{m}(\mathrm{h})}\left\{\mathrm{Z}({\mathrm{x}}_{\mathrm{i}}\right.)-{\left.\mathrm{Z}({\mathrm{x}}_{\mathrm{i}}+\mathrm{h})\right\}}^2$(2)

After the calculation of the semivariance, the variogram model was fitted, and the best-fit model was used to describe the relationship between semivariance and the lag. Frequently the semivariance increases with increasing lag but often reaches a maximum at some distance. The semivariance at this maximum value is called a sill. The lag when the sill is reached is called range, and this indicates the limit of spatial dependence. Sample distances further than the range are spatially independent. The semivariance at lag 0 is called the nugget and may reflect measurement error or spatial variation on a scale smaller than the observed grid size. The percent nugget was calculated to show the degree of spatial dependence, as equation 3 (Goovaerts 1999):

(3) $\mathrm{Percent}\mathrm{nugget}=\frac{\mathrm{nugget}}{\mathrm{sill}}\times 100$(3)

The percent nugget takes a value in the range of 0 to 100, and when <25% indicates a strong spatial dependence, 25–75% indicates a moderate spatial dependence, and >75% indicates a weak spatial dependence (Cambardella et al. 1994). Then, ordinary kriging was used to map each of the soil properties, and ArcMap 10.6.1 (Esri Inc. Berkeley, California) was used to display the maps.

3. Results and discussion

3.1. Soil profile characteristics and soil formation

There was a significant difference in the morphological features among the four soil profiles in different topographical positions (Figure 2). Especially, there were exceptional differences in the topsoil colors and A horizon thickness, likely affected by the differences in SOM content. Among the four soil profiles, profile 2, located at the lowest elevation (78.9 m), had the darkest color (10 YR 2/1 ~ 10 YR 1.7/1; Table 1) and the thickest A horizon reaching a depth of 55 cm (Figure 2; Table 1). Profile 1, located at a slightly higher elevation (79.2 m), had a similar to slightly lighter soil color (10 YR 2/1 ~ 10 YR 2/2) and much shallower A horizon thickness reaching a depth of 30 cm. Profile 3, located at a higher elevation (79.8 m), had even lighter soil colors (10 YR 2/3 ~ 10 YR 2/2) and the A horizon depth reaching a depth of 32 cm. Profile 4, located at the highest elevation (82.0 m), had the lightest soil color (10 YR 3/3) and the A horizon depth reaching a depth of 32/37 cm. The soil color reflected the total carbon (TC) content of the Ap1 horizon. The Ap1 horizon of profile 2 had the highest TC content of 106 g kg^‑1, followed by 78.9, 64.5, and 43.6 g kg^‑1 in profiles 1, 3, and 4, respectively (Table 2). This shows that among the four profiles, the TC content of the Ap1 horizon was higher as the field elevation decreased (Figure 4). The A horizon thickness was much higher in profile 2 compared with the rest of the soil profiles, likely affected by the concave topography unlike other profiles, which could accumulate eroded materials from surrounding areas (Kirkels, Cammeraat, and Kuhn 2014).

Figure 2. Photographs of the four soil profiles surveyed at different topographical positions.

Table 1. Description of morphological properties of the soil profiles.

	Depth (cm)	Moist soil color†	Mottling and concretion‡	Gravel§	Texture¶	Dominant parental material#
Profile 1
Ap1	0 - 12	10YR2/1	N	N	L	–
Ap2	12 - 30	10YR2/1	N	N	L	–
2 BA	30 - 40	10YR2/3	N	N	L	Ta-c
2B	40 - 49	7.5YR4/6	N	N	L	Ta-c
3B1	49 - 65	10YR4/6	N	N	L	En-a
3B2	65 - 81	10YR4/6	N	N	L	En-a
4B1	81 - 95	10YR5/6	N	N	CL	Paleosol
4B2	95 - 115+	10YR5/6	N	N	CL	Paleosol
Profile 2
Ap1	0 - 13	10YR2/1	N	N	L	–
Ap2	13 - 31	10YR2/1	N	N	L	–
2A	31 - 37	10YR1.7/1	N	N	L	Ta-c
2ACg	37 - 46	10YR3/2 (95%) 10YR4/3 (5%)	F TU-FeM (5%)	N	L	Ta-c
3ACg	46 - 55	10YR3/3 (95%) 10YR6/3 (5%)	F TU-FeM (5%)	N	L	Ta-d
4BAg	55 - 68	10YR4/4 (70%) 10YR5/6 (30%)	D TU-FeM (30%)	F SR S~LS (30%)	L	En-a
4Bg1	68 - 86	10YR6/6 (40%) 10YR4/6 (60%)	D TU-FeM (60%)	F SR S (30%)	CL	En-a
4Bg2	86 - 100+	10YR6/6 (60%) 10YR4/6 (40%)	D TU-FeM (40%)	F SR S (20%)	CL	En-a
Profile 3
Ap1	0 - 12	10YR2/3	N	N	L	–
Ap2/2A	12 - 32	10YR2/2	N	N	L	–
2 BA	32 - 42	10YR3/4	N	N	L	Ta-c
2B	42 - 52	10YR4/6	N	N	L	Ta-c
3B	52 - 68	10YR5/6	N	N	L	En-a
3Bg	68 - 90+	10YR5/6 (95%) 10YR4/6 (5%)	F TU-Fe (5%)	F SR S (40%)	L	En-a
Profile 4
Ap1	0 - 15	10YR3/3	N	N	L	–
Ap2	15 - 32/37	10YR3/3	N	N	L	–
2B	32/37 - 50	10YR4/6	N	N	L	Ta-c
3C1	50 - 66	10YR5/6	N	N	LS	En-a
3C2	66 - 100	10YR5/6	N	N	LS	En-a
4B	100 - 112+	10YR5/6	N	N	SL	Paleosol

†Based on the Munsell color chart.

‡N, none; F, faint; TU, tubular; FeM, Fe mottling; D, distinct.

§N, none; F, fresh; SR, subrounded; S, stone; LS, large stone.

¶L, loam; CL, clay loam; LS, loamy sand; SL, sandy loam.

#Ta-c, volcanic materials from Mt. Tarumae (2,500 ~ 3,000 years B.P.); Ta-d, volcanic materials from Mt. Tarumae (8,000 ~ 9,000 years B.P.); En-a, volcanic materials from Mt. Eniwa (19,000 ~ 21,000 years B.P.; Machida and Arai 2003).

Table 2. Results of physicochemical analysis of the soil profiles.

	Total carbon	Alo†	Feo	Sio	Alp	Fep	Allophane	Alo+½Feo	PAC
	(g kg^−1)							(%)	–
Profile 1
Ap1	78.9	26.5	9.94	9.14	5.89	3.75	77.0	3.15	1,780
Ap2	75.0	25.8	9.00	9.15	5.69	3.71	75.1	3.03	1,710
2 BA	43.4	76.4	21.5	34.9	5.07	.87	269	8.71	2,530
2B	15.9	56.9	15.0	28.9	3.03	.39	208	6.44	2,230
3B1	9.54	45.3	11.4	22.8	2.43	.33	165	5.10	2,000
3B2	6.94	38.4	9.22	21.0	2.11	.30	144	4.30	1,780
4B1	6.07	33.5	11.1	19.0	1.92	.40	127	3.91	1,810
4B2	6.51	34.6	10.7	20.7	1.48	.32	136	4.00	1,850
Profile 2
Ap1	106	24.2	8.59	6.62	9.35	4.67	55.4	2.85	1,940
Ap2	96.2	26.2	8.70	7.62	8.99	4.48	64.3	3.06	2,050
2A	105	24.9	8.51	6.96	9.16	4.66	58.5	2.91	1,970
2ACg	71.0	90.1	9.39	40.2	8.84	2.83	307	9.48	2,570
3ACg	44.7	79.2	10.9	37.1	5.91	.85	279	8.47	2,490
4BAg	29.1	34.4	10.8	14.3	4.32	1.64	113	3.97	2,080
4Bg1	14.4	29.9	9.68	13.6	3.08	.76	102	3.47	1,870
4Bg2	11.3	31.6	10.4	15.3	2.66	.54	111	3.68	1,810
Profile 3
Ap1	64.5	24.8	8.93	9.31	4.61	2.88	75.5	2.93	1,640
Ap2/2A	65.1	23.8	8.47	8.82	4.16	2.31	73.1	2.80	1,660
2 BA	39.3	46.9	14.3	19.9	4.77	2.53	158	5.41	2,280
2B	15.8	37.1	10.7	17.6	2.70	.54	131	4.25	1,940
3B	11.1	30.5	10.2	14.8	2.40	.47	108	3.56	1,760
3Bg	6.02	18.9	5.40	9.72	1.54	.34	68.1	2.16	1,150
Profile 4
Ap1	43.6	38.2	8.81	18.7	3.02	1.20	135	4.26	1,770
Ap2	45.1	39.0	9.29	19.5	3.09	1.21	139	4.37	1,790
2B	9.34	49.0	12.2	29.2	2.06	.28	192	5.51	1,910
3C1	5.22	35.6	10.7	22.3	1.66	.25	143	4.10	1,360
3C2	4.21	27.9	4.40	17.2	1.45	.27	111	3.01	1,220
4B	3.28	15.9	10.0	8.83	.80	.35	6.1	2.08	1,270

†Al_o, Al content extracted by acid oxalate; Fe_o, Fe content extracted by acid oxalate; Si_o, Si content extracted by acid oxalate; Al_p, Al content extracted by pyrophosphate; Fe_p, Fe content extracted by pyrophosphate; PAC, phosphate absorption coefficient.

The subsoils were very different among the four soil profiles, especially the presence of mottling, rock fragments, as well as the parent materials. Distinct tubular mottling was observed at 37 cm depth in profile 2, while faint tubular mottling was observed at 68 cm depth in profile 3, indicating seasonal high groundwater levels in these profiles, unlike in profiles 1 and 4, where no mottling was observed (Table 1). Subrounded rocks were also present in profiles 2 and 3, indicating that these rocks were deposited by water movement. These horizons were mixed with the parent materials of En-a, which were deposited 19,000 to 21,000 years B.P. from Mt. Eniwa (Machida and Arai 2003; Table 1). Therefore, the presence of watercourses is suggested in those eras in locations around profiles 2 and 3. Profile 1 was located at a lower elevation area compared with profile 3, but it lacked mottling and gravels. Therefore, the ancient water courses, as well as current groundwater movement, appear not to follow the present elevation in some areas of this field. Profile 4 was located at the highest elevation on top of the ancient sand dune formed by the En-a materials (Figures 1(b, c)). This location has likely experienced different moisture conditions compared with the other soil profiles. In general, the higher the moisture content is, the more the accumulation of SOM because of the lower microbial decomposition rate in oxygen-limited conditions (Saigusa, Shoji, and Otowa 1991).

All of the topsoil and subsoil horizons showed Al_o+½Fe_o values of more than 2%, and their morphological features were as follows: (i) Ta-c (volcanic materials from Mt. Tarumae, about 2,500 years B.P.) identifiable by its distinct soil texture and extremely high Al_o+½Fe_o content; (ii) Ta-d (volcanic materials from Mt. Tarumae, about 8,000 years B.P.) which has a distinct orange color after weathering and pumice-like texture; (iii) En-a (volcanic materials from Mt. Eniwa, about 18,000 years B.P.) identifiable by stratified soil texture from loamy to sandy due to water effects in the past; as well as some (iv) paleosols (Table 1; Machida and Arai 1992). The topsoil is assumed to be comprised of volcanic materials of more recent times, including Ta-a (volcanic materials from Mt. Tarumae, AD 1739), Ko-c₂ (volcanic materials from Mt. Komagatake, AD 1694), and Ta-b (volcanic materials from Mt. Tarumae, AD 1667), which were mixed by plowing in the Ap horizons. The Spfa-1 (volcanic materials from the Shikotsu caldera, about 40,000 to 45,000 years B.P.) was not identified within the 1-m depth in this field, though it is known to be deposited in the area (Kikuchi 2008). Horizons with Ta-c showed much higher allophane content compared with other horizons (Tables 1, 2). The morphological features of the four soil profiles evaluated in this study were similar to those reported on the middle terrace in this region (Kato et al. 2019; Kikuchi 2008; Tani et al. 2010).

The PAC values were also the highest, where the allophane content was high for the subsoils (Table 2). In the topsoil, PAC was high where allophane content was low such as the Ap1 and Ap2 horizons of profile 2; therefore, the contributions of several different soil components, such as the Al-humus complex on the PAC, need to be evaluated in more detail.

Soil texture was affected by the parent materials, and they were generally loam (Table 1). En-a is known to have distinctly different soil textures in locations where it was sorted according to the particle size by water after deposition (Machida and Arai 1992). In profiles 1, 2, and 3, En-a horizons were loamy to clay loam texture, but they were loamy sand in profile 4. This indicates that the coarse pumice-rich layer remained in profile 4, where it was not affected by water because of the higher elevation. The paleosols of profiles 1 and 4 were different in soil texture as well as allophane content, indicating that they came from different parent materials or had very different weathering processes due to the differences in the soil moisture conditions.

All of the soil profiles had PAC values of more than 1,500, which is the threshold used to classify Andosols in Japan according to the Soil Classification System of Japan (The Fifth Committee for Soil Classification and Nomenclature 2017). However, PAC values less than 1,500 were recorded in the 3Bg horizon in profile 3 and 3C1, 3C2, and 4B horizons in profile 4 (Table 2). Nevertheless, all of the soil horizons showed an Al_o+½Fe_o content of more than 2%, one of the requirements of andic soil properties in Soil Taxonomy (Table 2; Soil Survey Staff 2014) as well as the Soil Classification System of Japan. According to Soil Taxonomy, all of the profiles were classified into Andisols (udands), but there were two different types at the great group level. Profiles 1 and 2 were classified into Typic Melanudands because the cumulative thickness of the horizons, which showed that andic properties were more than 30 cm within a total of 40 cm thickness, with a color value of 2.5 or less and chroma of 2 or less. Profiles 3 and 4 were classified into Typic Hapludands. According to the Soil Classification System of Japan, all of the soils were classified into the Andosols great group because the thickness of horizons which showed Andic soil properties were 25 cm or more within 50 cm of the surface layer, but there were three types at the subgroup level. Profiles 1 and 3 were classified into Humic Allophanic Andosols because the topsoil had a color value of 3 or less and chroma of 3 or less, and the SOC content of the surface soils was 3% or more. Profile 2 was classified into Haplic Wet Andosols because the soil was assumed to be saturated in a certain period in a year, showing no gleyic characteristics, and oxidized Fe coating on the surface of peds and pores was observed. Profile 4 was classified into Brown-humic Allophanic Andosols because the surface layer had a SOC content of 3% or more, but the color did not qualify for humic characteristics.

3.2. Characteristics of the surface soil properties

The surface soil properties showed distinct characteristics of Andosols and presented exceptional variations. The mean TC content was high at 62.5 g kg⁻¹, with a huge variation ranging from 32.5 to 124 g kg⁻¹, with a CV of 35.2% (Table 3). This was more variable compared with the past studies conducted under Andosols. In the same study area, Niwa et al. (2019) reported TC contents varying from 29.6 to 82.4 g kg⁻¹ with a CV of 24.7% within a single field, where the parent materials were volcanic material overlying late Pleistocene deposits. The high variation of TC content in the surface soil found in our study might have resulted from the soil formation process undergone in this region; therefore, similar variations are expected in other fields in the Tokachi region. Kinoshita et al. (2016) studied the within-field variation of TC in a 100-ha coffee agroforestry farm with an elevation difference of 269 m in Costa Rica, where TC ranged from 48.1 to 172 g kg⁻¹ with a CV of 33.2%. These authors reported that andic soil properties and topography largely controlled the TC variation, which follows what was found in our study (discussed in detail in section 3.4).

Table 3. Descriptive statistics of selected surface inherent soil properties.

	Total carbon	Alo†	Feo	Sio	Alp	Fep	Allophane	Alo+½Feo	CEC	PAC
	g kg^−1							%	cmolc kg^−1
Mean	61.5	33.4	12.1	14.5	5.88	2.97	107	3.94	30.5	1,830
Minimum	32.5	20.0	7.72	5.28	2.76	.70	39.0	2.39	20.3	1,310
Maximum	124	48.8	20.1	25.0	12.9	9.87	214	5.67	55.7	2,520
SD‡	21.6	5.82	1.97	4.41	2.28	1.70	29.3	0.64	7.01	178
CV§ (%)	35.2	17.4	16.3	30.4	38.8	57.3	27.4	16.2	23.0	9.70

†Al_o, Al content extracted by acid oxalate; Fe_o, Fe content extracted by acid oxalate; Si_o, Si content extracted by acid oxalate; Al_p, Al content extracted by pyrophosphate; Fe_p, Fe content extracted by pyrophosphate; CEC, cation exchange capacity; PAC, phosphate absorption coefficient.

‡SD, standard deviation.

§CV, coefficient of variation.

Al_o represents the Al contained in both noncrystalline and paracrystalline clay minerals (i.e., allophane, imogolite) as well as the Al-humus complex (Blakemore, Seatle, and Daly 1987). The Al contained in the Al-humus complex is extracted as Al_p. The variation of Al_o was much less, with a CV of 17.4%, compared to Al_p, with 38.8%. There was a significant positive relationship between TC and Al_p (r = 0.94; Table 4), which suggests that Al-humus complex is one of the key factors for the accumulation of SOM. Allophane content varied from 39.0 to 214 g kg⁻¹ with a CV of 27.4%, slightly less variable compared with TC (Table 3), indicating that clay mineralogy varies less compared to TC in this field. There was a significant negative correlation between allophane and Al_p as well as TC (r = −0.69; Table 4). This is because the formation of allophane and imogolite is limited by aqueous Al activity, and the formation of the Al-humus complex is known to lower the aqueous activity of Al (Dahlgren, Shoji, and Nanzyo 1993). Also, soluble humic substances are known to perturb the interaction between hydroxyl-Al ions and orthosilicic acid, which hinders the formation of imogolite and allophane (Inoue and Huang 1990).

Table 4. Pearson correlation coefficient of selected surface inherent soil properties (p < 0.001, except that indicated accordingly).

	Total carbon	Alo	Feo	Sio	Alp	Fep	Allophane	CEC
Alo†	−0.54
Feo	−0.55	0.55
Sio	−0.81	0.82	0.66
Alp	0.94	−0.56	−0.51	−0.82
Fep	0.88	−0.50	−0.53	−0.76	0.85
Allophane	−0.69	0.94	0.60	0.90	−0.72	−0.64
CEC	0.94	−0.53	−0.49	−0.77	0.91	0.82	−0.68
PAC	0.66	n.s. ‡	−0.16	−0.22	0.59	0.60	−0.14	0.66

†Al_o, Al content extracted by acid oxalate; Fe_o, Fe content extracted by acid oxalate; Si_o, Si content extracted by acid oxalate; Al_p, Al content extracted by pyrophosphate; Fe_p, Fe content extracted by pyrophosphate; CEC, cation exchange capacity; PAC, phosphate absorption coefficient.

‡n.s, not statistically significant.

CEC had a mean of 30.5 cmol_c kg⁻¹ and ranged from 20.3 to 55.7 cmol_c kg⁻¹ with a CV of 23.0% (Table 3). It is known to be affected by the type and amount of both clay minerals and humic substances (Akbas, Gunal, and Acir 2017; Manrique, Jones, and Dyke 1991), and TC was the most correlated property in our study (r = 0.94; Table 4), indicating that the amount of humic substances controls the variation of CEC. This also means that many of the negative charge sites are pH-dependent in this field (Nanzyo, Dahlgren, and Shoji 1993). There have been limited studies done to analyze the within-field variation of CEC under Andosols. The CEC of this study was extremely high compared with past studies in upland fields of various soil types. In nine fields across four states in the U.S.A., the mean CEC of each field ranged from 2.80 to 24.0 cmol_c kg^−1, and the maximum values ranged from 3.50 to 32.9 cmol_c kg⁻¹ (Kweon, Lund, and Maxton 2013). In their study, the CV was lower than our study except for just one field, which had a CV of 32.6%, indicating that our field showed a very high variation for an upland field. CEC is an extremely important indicator of soil’s capacity to retain base cations such as Ca²⁺, Mg²⁺, and K⁺ important for crop productivity. In Hokkaido, the optimum base saturation is recommended to be between 60% and 80% (Department of Agriculture, Hokkaido Government 2020). When the CEC values have variation, the base saturation is expected to be also variable, and it can be below the optimum range, especially in areas with high CEC, or can be excessive in areas with low CEC. These can potentially lower crop productivity when the supply of the base cations is limited. The CEC was strongly correlated with TC (r = 0.94), indicating that most of the negative charge sites are pH-dependent, which are known to have specific adsorption for trace metals (McBride 1994; Mouni et al. 2009, 2013; Nanzyo, Dahlgren, and Shoji 1993). Therefore, in high CEC areas, the essential plant nutrients such as Cu²⁺ and Zn²⁺ may become limited, while potential pollutants such as Pb²⁺ and Cd²⁺ can be retained.

PAC had a mean value of 1,831 and ranged from 1,310 to 2,520, with a CV of 9.70% (Table 3). Seven samples out of 530 samples had PAC values less than 1,500, the value used to classify Andosols by the Soil Classification System of Japan (The Fifth Committee for Soil Classification and Nomenclature of Japanese Society of Pedology 2017). However, this value can be replaced with an Al_o+ ½Fe_o value of more than 2%, and all the 530 samples satisfied the criteria. Since the soil samples were taken from 0 to 25 cm depth, the entire field can be classified into Andosols. PAC is a measure of soil’s capacity to fix P through mechanisms such as adsorption or precipitation, and higher values indicate higher fixation capacity. The variation was much less compared with TC or allophane in this field, and most of the PAC values ranged from 1,500 to 2,000. PAC is affected by non- and paracrystalline clay minerals (i.e., allophane and imogolite), Al-humus complex, ferrihydrite, Al hydroxide, exchangeable Al, Ca, and Mg (Editorial Committee for Methods of Soil Environmental Analysis 1997). The PAC values were positively correlated with TC (r = 0.66; Table 4) and negatively correlated with allophane (r = ‑0.14), showing that the Al-humus complex is likely the dominant factor influencing PAC in this study. Niwa et al. (2019) also found a positive correlation between TC content and PAC (R² = 0.70) in the 40 fields surveyed in this region, where the parent material was a mixture of volcanic material and alluvial deposits. However, Hashimoto et al. (2012) showed Al_o to be the most important predictor of PAC for allophanic Andisols when they studied 1,300 soil samples collected from nationwide Japan. Our field was classified as allophanic Andosols, but the correlation between PAC and Al_o was low and not significant (Table 4), which contradicted their finding. This might have been caused by the fact that in our study the Al_o range was smaller (20 to 48.8 g kg⁻¹) compared to that found by these authors (1.2 to 99.2 g kg⁻¹). Also, in our study, the PAC was strongly correlated with CEC (r = 0.66) while no correlation between these properties was found by Hashimoto et al. (2012). If most of the negative sites on the surface of the clay colloids are occupied by base cations such as Ca²⁺ and Mg²⁺, it is expected that these might control the PAC in the field (Nanzyo, Uhm, and Shoji 1992). Although the relation between PAC and base cations was not evaluated in this study, we hypothesize that this could be the factor controlling the PAC in this field. This implies that the factors affecting PAC are more complex at a within-field scale and require further studies since PAC is very important when considering optimum phosphate fertilizer management.

3.3. Spatial structure of the surface soil properties

The spatial structure of each measured property was assessed using geostatistics. A spherical model was fit for all properties (Table 5) and then used to map the soil properties in ordinary kriging (Figure 3). In this study, most of the soil properties showed strong spatial class except for Al_o, Fe_o, and PAC (Table 5), showing clear zones in the maps (Figure 3). Strong spatial class indicates that the soil properties are controlled by intrinsic variations of soil characteristics (i.e., texture and mineralogy) rather than extrinsic variations (i.e., fertilizer and tillage; Cambardella et al. 1994). All of the soil properties measured in this study are regarded as inherent soil properties that are intrinsic properties of soils. The strongest spatial class was identified for TC, CEC, and Al_p with 0% nugget. These properties had a long range with 226, 213, and 188 m, respectively, indicating that the values change gradually over a long distance and presenting large zones (Figures 3(a, b)). These properties show high correlation coefficients with each other (Table 4), and they also show similar spatial structures. The high values were concentrated in the middle of the field where the elevation was lower (Figures 1b, 3(a, b)), and three distinct zones of lower, higher, and lower values were apparent from the north to the south following the microtopography formed due to the two ancient sand dunes existing in this field. Kato et al. (2019) reported a moderate spatial dependence of TC (% nugget of 44), although the TC variation was higher compared to our study, which contradicts our findings. The strong spatial dependence of TC in our study can be explained by the fact that the sampling interval was relatively smaller (12-m sampling interval, equivalent to 80 samples ha⁻¹) compared to 5.8 samples ha⁻¹, used by Kato et al. (2019). The TC nugget in our study was zero compared to 449 reported by the same authors, showing that both the sampling interval and the TC measurement were appropriate in our study (Goovaerts 1998). TC values are also known to be correlated with hot water extractable N, which can be used as an indicator of available N in the soil (Curtin et al. 2006; Niwa et al. 2016). The spatial distribution of soil available N in this field is expected to be similar to the TC. Thus, its spatial management could be accomplished using the TC map. This also suggests that TC spatial variability maps can be used to define N management zones at the regional scale. Fe_p was also correlated with TC (r = 0.88), and it also showed a relatively long range at 170 m (Table 5). Allophane, Al_o, and Si_o showed relatively shorter ranges between 85.5 and 122 m (Table 5). The values were low in the center of the field where the elevation was lower, and TC was high (Figures 1b, 3(c,d)). These were expected because of the negative relationship between TC and allophane discussed previously in section 3.1. However, there was a more complex spatial variation for allophane compared with TC, especially in the southern part of the field around the ancient sand dune.

Figure 3. Kriged maps of selected inherent soil properties on the surface soil.

Table 5. Parameters of the semivariograms for selected surface inherent soil properties.

	Nugget (C0)	Sill (C0+C)	Nugget %‡	Range (m)	Spatial Class§
Total carbon	0.00	710	0.00	226	Strong
Alo†	8.46	29.0	29.2	85.5	Moderate
Feo	1.87	2.75	68.0	168	Moderate
Sio	0.84	22.6	3.72	122	Strong
Alp	0.00	7.30	0.00	188	Strong
Fep	0.43	3.33	12.9	170	Strong
Allophane	133	898	14.8	118	Strong
CEC	0.00	67.8	0.00	213	Strong
PAC	8,534	24,778	34.4	139	Moderate

†Al_o, Al content extracted by acid oxalate; Fe_o, Fe content extracted by acid oxalate; Si_o, Si content extracted by acid oxalate; Al_p, Al content extracted by pyrophosphate; Fe_p, Fe content extracted by pyrophosphate; CEC, cation exchange capacity; PAC, phosphate absorption coefficient.

‡Nugget % = C₀/(C₀+C)×100.

§Strong, Nugget % ≤ 25; Moderate, 25 < Nugget % ≤ 75.

PAC showed a moderate spatial structure with a relatively shorter range at 139 m (Table 5). Although it was most correlated with TC, the spatial structure and the map were different (Figures 3(a,e); Table 5). The high TC area and high PAC area were similar in the center of the field; however, there were areas of high PAC and low TC that existed on the southwestern side of the field. This information warrants further studies to understand the complex factors affecting PAC and how it varies spatially within a field.

3.4. Factors affecting the spatial variation of the surface soil properties

The spatial variation of TC was very high in this field (CV = 35.2%; Table 3), and it also had a strong spatial class with zero nugget percent (Table 5). There was a strong negative relationship between elevation and TC (r = −0.73; Figure 4). The change in TC content was much higher from the elevation of about 79.0 to 80.0 m, while much less TC change was observed where the elevation was more than 80.0 m. The TC content and elevation data of the Ap1 horizon of the four soil profiles were also plotted (Figure 4). A hypothesis can be formulated that the lower elevation site has a higher groundwater level that limits the availability of oxygen in the soil profile and, in turn, reduces the mineralization of SOM (Saigusa, Shoji, and Otowa 1991). This hypothesis is true if profile 1 is ignored and when comparing only profiles 2, 3, and 4. Profile 2 was located at the lowest elevation (78.9 m) and had the shallowest mottling at 37 cm depth, while profile 4 was at the highest elevation (82.0 m) and lacked mottling in the entire soil profile. If elevation has some associations with the depth of the seasonal water table in this field, the elevation between 79.0 and 80.0 m determines whether shallow mottling within 50 cm of a soil profile is present or a complete lack of mottling. The significance of this impact may appear in the large change of TC in this elevation range (Figure 4). Above 80.0 m, the change in the depth of the seasonal water table may be deeper than 1 m in the soil profile, and this difference will not affect much of the SOM decomposition. However, this hypothesis does not apply to profile 1, where the TC content was the second highest but lacked mottling within 1-m of the soil profile (Figure 1; Table 1). The Ap1 horizon of profile 1 had the highest Al_o content among profiles 1, 2, and 3 and also had very high allophane content. In this profile, Al that could react with humic substances was more available, especially compared with profile 3. This Al could form an Al-humus complex that can reduce microbial decomposition of complexed humic acid (Saigusa, Shoji, and Otowa 1991). Therefore, the combination of moisture content and the mineralogical composition of the soil profile are thought to be the reason behind the spatial variation of TC, which was overall influenced by microtopography. Our findings highlight the role of topography in describing the spatial variability of inherent soil properties in this region.

Figure 4. The relationship between elevation and total carbon for both surface as well as the Ap1 horizon of the soil profiles.

The spatial variation of the surface allophane content was negatively correlated with TC content (Table 4); however, the profile characteristics should also be considered for further evaluation. Profile 1 had a higher allophane content in the Ap1 horizon compared with profile 3, even though TC content was higher in profile 1 (Table 2). The surface allophane content appears to be affected by humic substances but also by plowing and mixing of the allophane-rich subsurface layer into the topsoil. Therefore, the spatial dependence of allophane was much less compared with TC and had a shorter range of 118 m (Table 5). These complex effects of intrinsic mechanisms of inhibition of allophane synthesis by humic substances and anthropogenic influence of plowing caused the spatial dependence to be weaker.

Overall, this field had two ancient sand dunes that formed microtopography. This microtopography affected the course of ancient rivers and also changed the depth of the current seasonal water table. The microtopography effect, combined with the presence of the depth of an allophane-rich subsurface horizon, appears to dominate the spatial variation of the surface’s inherent soil properties. The spatial variation of TC dominates the variation of CEC, while the variation of PAC appears to be affected by factors other than both TC and allophane; thus, further investigations are needed.

4. Conclusion

Exceptional spatial variations of inherent soil properties were found in an upland field in the Tokachi district, where TC varied more compared with clay minerals. The microtopography affected the within-field spatial distribution of soil properties by altering the depth of the seasonal water table, which is a major factor affecting the TC content in this region. Volcanic materials of different ages were found in the soil profiles, and their distributions were also dependent on the topographical position. Allophane content in the soil surface showed a weaker spatial dependence compared with TC, which was affected by both the TC content but also anthropogenic factors such as plowing by incorporating the allophane-rich subsurface layers into the topsoil. Agronomically important properties such as CEC were strongly affected by the TC content, while PAC showed weaker correlations to both the TC content and clay minerals.

To optimize nutrient management in the large-scale fields in the Tokachi region, management zones need to be created, taking into account the within-field microtopography. Base cations can be managed based on the CEC map which had a strong correlation to the TC map. The management of N can also be based on the TC map. The estimation of the TC content can be conducted by using remote sensing equipment such as an unmanned aerial vehicle or a satellite, and the estimation model can also incorporate the elevation data. The estimation of PAC is more challenging since both TC and clay minerals affect the property, and the relationship to the elevation was weaker. A future study needs to be conducted to assess the relationship between these inherent soil properties and nutrient availability. Also, the mapping of these soil properties using remote sensing should be attempted to alleviate the cost and time required to obtain the soil maps.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was funded by Loginet Japan Co., Ltd., JSPS Kakenhi Grant number JP22K05590, and Obihiro University of Agriculture and Veterinary Medicine. It also had technical support from Dr Shunzo Kawajiri from the Kyushu Institute of Technology and Dr Atsuhiro Saegusa from the Kitami Institute of Technology.