Maria Eliza TUREK，Robson AndréARMINDOand Ole WENDROTH

1 Graduate Program of Environmental Engineering(PPGEA)，Federal University of Paraná(UFPR)，Curitiba 81531 980(Brazil)

2 Depar tment of Physics(DFI)，Federal University of Lavras(UFL A)，Mailbox 3037，Lavras 37200 900(Brazil)

3 Depar tment of Plant&Soil Sciences，University of Kentucky(UK)，Lexington 40546-0091(USA)

ABSTRACT Soil physical quality(SPQ)assessment is an important part in the evaluation of soil use，management，and conservation.It can be assessed using several physical properties，hydraulic indices，and functions.Soils from tropical and temperate regions represent different physical behaviors，and the quantif ication of their physical properties is important to support soil evaluation and modelling.The objective of this study was to evaluate the SPQ in a subtropical f ield under maize crop cultivation according to its physical properties，hydraulic indices，and functions in an attempt to infer the spatial variability and to determine the behavior of soil physical structure across a spatial domain.Commonly used soil key physical variables，such as texture，bulk density，total porosity，saturated hydraulic conductivity，and organic carbon content，were measured in a regular grid with a soil sampling density of 30 points per hectare，covering an area of 0.5 ha.Saturated hydraulic conductivity varied strongly between subsamples and in the f ield，suggesting the heterogeneity of the soil structure regarding water drainage.The physical variables were combined with other indicators，which were based on the soil water retention curve and the pore size distribution(PSD)function.Correlation analysis was performed to verify the relationship between the measured and calculated variables，and some strong linear correlations were revealed，such as between aeration energy index and microporosity(r=0.608)and water retention energy index with microporosity(r=0.532)and with bulk density(r=0.541).For most sampled locations，the shape and location parameters of PSD showed results outside of the optimum ranges，whereas the hydraulic energy indices and cumulative hydraulic energy functions presented values that were similar to those found for some tropical soils described in the literature.The spatial variability of these indices was described using semivariograms and kriged maps，indicating the variability of the SPQ in this f ield.

Key Words: geostatistics，integral energy，pore size distribution，soil physical indicators，soil physical quality，soil water retention curve

INTRODUCTION

Soil physical quality(SPQ)is related to the ability of the soil to promote satisfactory conditions for crop growth and development and to buffer the effects of extreme weather conditions.In several studies，SPQ indicators that reveal the interaction among water inf iltration and storage，heat，gas exchange，and organic carbon dynamics have been proposed(Shahabet al.，2013;Armindo and Wendroth，2016;Reynolds，2017).These indicators are tools to reveal variations and patterns in time and space with substantial efforts spent in the previous years to identify threshold values for soil physico-chemical properties，such as organic carbon content，dry bulk density，and hydraulic variables(Reynoldset al.，2009).

Soil structure directly inf luences SPQ because it regulates processes，such as water inf iltration，redistribution，and retention，gas exchange，soil organic matter，nutrient dynamics，root penetration，and susceptibility to erosion(Rabotet al.，2018).Therefore，many efforts have been made to evaluate soil structure，usually based on some key physical variables，such as texture，bulk density(ρb)，macroporosity(macro)，microporosity(micro)，and total porosity(φ)，saturated hydraulic conductivity(Ks)，and organic carbon(OC)content.These studies revealed soil structural dynamics over time and space，such as the effect of soil compaction，reduced tillage(Schlteret al.，2018)，mulching schemes，irrigation，drainage，water diversion structures for soil conservation(Jourgholamiet al.，2018)，and natural weather conditions(Bryket al.，2017).The easy-to-measure properties，ρbandφ，are indicators of soil structure and are contained in most soil databases to provide information for preliminary evaluation of SPQ(Hkansson and Lipiec，2000;Rabotetal.，2018).Saturated hydraulic conductivity is another commonly measured property that carries hydraulic information ref lecting the interaction between soil structure and biological activity，root growth，and swelling and shrinking dynamics(Perkinset al.，2007;Sandinet al.，2017).

Besides those key properties，two widely used hydraulic functions that carry essential information are the soil water retention curve(SWRC)(Ghanbarianet al.，2015)and its f irst derivative，the pore size distribution(PSD)function(Reynoldset al.，2009).Because these functions and the geometry of pores are strongly related to SPQ，changes in soil structure are manifested by differences in both PSD and SWRC(Levy and Mamedov，2002).One limitation to their wide-spread use is that their measurement is time-consuming and costly，especially for large-scale descriptions(Van Looyet al.，2017).

Although measured SWRCand PSDover the entire range of soil water content is rather valuable to the characterization of soil physical structure，most SWRC-and PSD-based approaches do not use entire curves to evaluate SPQ，but instead use only a few points.In contrast，other studies suggest that SPQ assessment based on the entire SWRC and PSD is more appropriate(Minasny and McBratney，2003;Shiet al.，2011;Armindo and Wendroth，2016).According to Armindo and Wendroth(2019)，it was most likely Mualem(1978)who f irst def ined the energy—work required to drain a volume of saturated soil by integrating the SWRC between inf inity and zero limits of soil water pressure head.Later，Minasny and McBratney(2003)integrated the Van Genuchten-Mualem(VGM)model between 15 000 and 100 hPa to calculate the integral energy，per unit of plant-available water(PAW)，held by the soil in this range.The concept of the SWRC integral was also applied with distinct SWRC models，indices，and approaches to characterize SPQ(Asgarzadehet al.，2011;Shietal.，2011;Chamindu Deepagodaetal.，2013;Armindo and Wendroth，2019).

New experimental designs and several pedotransfer functions(PTFs)have been developed to estimate SWRCparameters with less cost than the traditional laboratory measurements to associate hydraulic functions with easier-to-measure properties(Nemeset al.，2006;Huf dos Reiset al.，2018).The applied method should be carefully chosen because statistical PTFs depend on the calibration dataset;some are based only on texture or non-structural properties，implying substantial uncertainty(Castellini and Iovino，2019).For example，the Beerkan method for estimating soil transfer parameters(BEST)(Lassabatreet al.，2006)includes a f ield inf iltration experiment where the shape parameters of the SWRC and PSD function are based on soil textural information.A different approach was presented by Prevedello and Armindo(2016)，wherein the inf iltration was estimated based on the SWRC parameters andKs.The BEST method was recently applied by Castelliniet al.(2018)to estimate PSD and evaluate SPQ in the f ield under winter wheat and by Di Primaet al.(2018)for soils cultivated with citrus orchards under several tillage systems.Soil structure has been evaluated by means of PTFs，as shown by Fentonet al.(2017)and Huf dos Reiset al.(2018)，who applied Rosetta(Schaapet al.，2001)and Splintex(da Silvaet al.，2017b)to obtain PTF-derived SPQ indicators，as well as by da Silvaet al.(2017a)，who evaluated the spatial variability of soil key physical variables based on the results from Splintex and Rosetta.Bacheret al.(2019)tested the inf luence of different schemes for obtaining the SWRC on SPQ，including the use of Rosetta，f inding good estimates of PAW.However，as for any method of estimation，uncertainties are associated with the BEST method，as reported by Castelliniet al.(2018)，as well as with Rosetta and Splintex(da Silvaet al.，2017a;Huf dos Reiset al.，2018).

Armindo and Wendroth(2016)examined f ive hydraulic energy indices and two cumulative energy functions，based on the integration of the SWRC，to evaluate SPQof soils from different regions of temperate climates，such as Harzgerode and Kraatz in Germany and Lexington in the USA，and tropical climates，located near the Brazilian cities of Boca da Mata-AL，Maringá-PR，Araponga-MG，Choró-CE，and Alegre-ES.The results revealed differences in SPQ because of the type，use，and management of these soils，manifesting their structural behavior.The authors suggested that the energy indices and functions should be studied in other soils to establish threshold-ranges or at least provide an estimation of the magnitude of their values in some target applications，e.g.，agriculture.The objective of this work was to evaluate SPQ indices and their spatial dependence in a clayey Cambisol under maize crop.

MATERIALS AND METHODS

Studied site and soil sampling

The f ield site has soil classif ied as Cambisol and was located in the Applied Agricultural Sciences Farm of the Federal University of Paraná(UFPR)，near Pinhais City，Paraná，Brazil(25°23′20′′S，49°7′50′′W).Disturbed and undisturbed soil samples were collected at 15 sampling points from the 0—0.1 m soil layer using a regular 25 m×25 m grid，covering an area of 0.5 ha.This soil layer was chosen because of its effect on many environmental and agronomic processes，such as seed germination，crop growth，soil erosion，runoff，and water inf iltration(Cherobimet al.，2015).At each sampling point，three undisturbed subsamples were taken as replicates，resulting in a total of 45 samples.The undisturbed soil cores were collected using thin-walled stainless steel cylinders(3 cm in diameter and 4 cm in height)and wrapped in a plastic f ilm until laboratory analysis.

Laborator y analysis

The disturbed soil material was dried(60°C for 48 h)and sieved through a 2-mm mesh for texture and OC content determination.The textural analysis was conducted based on the Bouyoucos hydrometer methodology(Gee and Or，2002)and the soil type was then established according to the United States Department of Agriculture textural triangle(Soil Survey Staff，2014).The OC content was measured by cold digestion of Na2Cr2O7(Sampaioet al.，2012).

The intact soil cores were prepared for analysis using serrated knifes and paint brushes to prevent sealing of the surface，which might have otherwise been created if smooth blades were used.After the initial preparation，the core samples saturated by capillary rise set the pressure boundary condition at the bottom by establishing a water table at 0.02 m.Saturated soil water content(θs)was then measured using an electronic balance(BEL model S3102，Bel Engineering，Italy)with a precision of 0.01 g.

SoilKswas measured in a constant head permeameter using a Mariotte tube to control the pressure head over the saturated soil cores，similar to the one described by Reynolds(2009).For each soil core，three replicate measures of water f low were taken after reaching a pressure equilibrium with the Mariotte tube to obtain meanKsbased on nine values.Soil cores were then placed on a pressure plate apparatus for determination of the volumetric soil water content(θ)at 0.6-m pressure(θ(0.6 m))and afterwards dried in an oven at 105°C for 48 h.The measuredθ(0.6 m)was set equal tomicrobecause，at this pressure head，all pores with equivalent diameter larger than 50μm were drained.Macroporosity was calculated by the difference between measuredθsandmicro;θ(0.6 m)was also chosen as an input parameter for the Splintex model to build SWRCs，following da Silvaet al.(2017b)，who found the best SWRC results usingθ(0.6 m)after testingθ(0.3 m)，θ(0.6 m)，andθ(1.0 m).Bulk density was calculated as the ratio between the oven-dry(105°C)soil mass and the cylinder’s volume.Solid particle density(ρs)was determined by measuring the volume of 20 g oven-dried and sieved soil material in alcohol.

SWRC

The SWRC was characterized by the relationship betweenθand soil water pressure head(h)and could be represented by the VGM model:

whereθ(h)is the volumetric soil water content ath;θris the residual soil water content;andα，n，andmare f itting parameters，withm=1-1/nandh＞0 for unsaturated conditions.

The model Splintex was applied to every sampling location to estimate the parametersαandnof Eq.1(da Silvaet al.，2017b).Splintex is a physically based model composed of two PTFs，with the advantage that in case of lack of calibration，it can be directly based on results obtained for a wide range of soils.The input variables required for Splintex are other soil physical properties，such as texture，ρb，ρs，and optionally，θs，and any other measured point on the SWRC.For more details，see da Silvaet al.(2017a，b)and Huf dos Reiset al.(2018).

Preliminar y experiments designed to build the SWRC

To derive a procedure for estimating the SWRC，we designed two preliminary simulations using a clayey and a silty clay soil randomly chosen from the Unsaturated Soil Hydraulic Database(UNSODA)(Nemeset al.，2001)，both considered being close in soil texture to the ones evaluated in this study.Splintex was then applied to obtain the SWRC based on the particle size distribution，ρb，ρs，φ，andθ(0.5 m).We usedθ(0.5 m)as input information because it was the closest measured value toθ(0.6 m)，which would be further used in this work.

The f irst simulation，hereafter referred to as Approach A，was performed by setting the parameterθstoφand using the outputsθr，α，andnprovided by Splintex as recommended by da Silvaetal.(2017a).Then，in a second simulation，hereafter referred to as Approach B，the largest measured value ofhin the SWRCs was chosen and assumed to correspond toθr.In both soils，the larger measured point wasθ(158.5 m);therefore，we considerθr≈θ(158.5 m).Hence，for Approach B，we used only the values of parametersαandnprovided by Splintex because the measured values were set asθs=φandθr=θ(158.5 m).

Obtaining the SWRCs

The value ofθ(158.5 m)was measured using a dewpoint meter(Decagon Devices，USA)，which is a widely used approach for measuringhat low soil water contents.Using the disturbed soil samples，hwas measured several times until a mean value was obtained to represent the entire f ield.The measured gravimetric soil water content(θg)at 158.5-m pressure(θg(158.5 m))for these samples was multiplied byρb/ρwater，whereρwateris the density of water，of each sampling point(Table I)to obtain the volumetricθ(158.5 m).Splintex was then fed with the data from the cumulative particle size distribution，described for particle diameters smaller than 0.002 mm(clay)，between 0.002—0.05 mm(silt)，and between 0.05—2.0 mm(sand)，ρb，ρs，and the measured pointsθsandmicroatθ(0.6 m)of the SWRC(Table I).

Soil physical indicators and functions

The indices of air capacity(AC)and PAW，derived from SWRC，and the median，mode，and mean of pore diameter(dmedian，dmode，anddmean，respectively)，kurtosis，skewness，and standard deviation(SD)，derived from PSD(Reynoldset al.，2009)，were calculated based on the parameters of Eq.1.The index ACwas calculated as AC=θs-θ(1.0 m)and PAW as PAW=θ(1.0 m)-θpwp，whereθpwpis the permanent wilting point，obtained ath=150 m.Field capacity(θfc)was obtained from the dynamic criterion(Prevedello，1999;Tureketal.，2019)，which considers that this point is reached when the soil drainage f lux equals to a small percentage(1%)ofKs(Armindo and Wendroth，2016;Tureket al.，2019).Then，the hydraulic energy indices of aeration absolute energy(Aa)，absolute water retention energy(WRa)，relative aeration energy(Ar)，relative water retention energy(WRr)，and relative air-water energy(AWr)and hydraulic energy cumulative functions examined by Armindo and Wendroth(2016)were investigated.The absolute energy indicesAaand WRaare results of the numerical integration of the VGM equation:

TABLE I Soil texture(gravimetric content)，bulk density(ρb)，saturated hydraulic conductivity(K s)，the standard deviation(SD)of K s，organic carbon(OC)content，microporosity(micro)，and macroporosity(macro)of 15 sampling points in a subtropical soil under maize crop cultivation in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil

The relative indices(Ar，WRr，and AWr)were calculated as:Ar=Aa/(Aa+WRa)×100，WRr=WRa/(Aa+WRa)×100，and AWr=Aa/WRa×100.The cumulative functions were calculated by the integral of the SWRC with a gradual increase ofθ.

Statistical and geostatistical analyses

The measured variables were evaluated according to mean values，SD，and coefficient of variation(CV)at each point and summarized for the entire area.The soil propertyKsf its better to a log-normal than a normal distribution.Therefore，mean values were calculated as the geometric mean.

The mean SWRC curve for the entire area was calculated using the mean values of the VGM parameters.The difference between this curve and those obtained at each of the 15 locations was evaluated by the root mean square error(RMSE)，comparing values ofθbetweenθ(0)andθ(150 m)(Morenoet al.，2019).

The geoR package(Ribeiro and Diggle，2018)was used to generate the spatial distribution of the analyzed key variables(Grego and Vieira，2005;Gülseret al.，2016)applying the Box-Cox parameter(λ)transformation when necessary.Experimental semivariogram data were f itted to Gaussian，spherical，Maternκ=2.5，and exponential models to choose the most parsimonious model by means of the maximum loglikelihood(VER).The selected model was used to krige maps and provide the parameters of practical range(a)，nugget(C0)，and partial sill(C1)(Kamimuraet al.，2013).The degree of spatial dependence(DSD)(DSD=C0/(C0+C1)×100)was applied to the variables associated with SPQ，according to the criterion presented by Cambardellaet al.(1994).

RESULTS AND DISCUSSION

SWRC

The results of approaches A and B are presented in Fig.1，together with the measured(UNSODA)data(Nemeset al.，2001).The insertion of a large value of soil water pressure head to representθrimproved the performance of Splintex in representing the measured SWRC.For the clayey soil，with the same textural composition as the studied Cambisol(Table I)，RMSE between measured and estimated curves was 0.101 8 m3m-3(Approach A)，and it decreased to 0.012 9 m3m-3when Approach B was applied.For the silty clay soil，the measures of RMSE were 0.099 6 m3m-3(Approach A)and 0.039 5 m3m-3(Approach B).The insertion of a measured point asθrtogether with measuredθsimproved the representativeness of the SWRCs both visually and statistically.Other studies have shown that the insertion of measured points has greatly reduced the deviations between the estimated VGM-parameter-based and the measured SWRCs(da Silvaetal.，2017a，b;Huf dos Reiset al.，2018).

Fig.1 Soil water retention curves obtained using a clay soil(a)and a silty clay soil(b)from the Unsaturated Soil Hydraulic Database(UNSODA)(Nemes et al.，2001)via approaches A and B.In Approach A，saturated volumetric soil water content(θs)was set to total porosity(φ)，and residual volumetric soil water content(θr)and the f itting parametersαand n of the Van Genuchten-Mualem model from Splintex were used.In Approach B，θs=φ，θr was set to volumetric soil water content at 158.5-m pressure(θ(158.5 m))，andαand n from Splintex were used.

Considering the previous simulations，the f inal parameter conf iguration for Eq.1 in the studied Cambisol wasθs(measured saturated water content)，θr(measuredθ(158.5 m))，andαandnas output from the model Splintex(Table II).

The individual graphical mean representation of each sampling point is shown in Fig.2，where the one that behaved closest to the mean curve(RMSE=0.011 5 m3m-3)of the studied f ield is highlighted as a bold solid line(Point 14).The estimated f ield variability of the SWRC across the 15 points in the studied area is presented in Fig.2.The natural spatial variability associated with inevitable experimental and methodological errors，including measurements taken in the laboratory，sampling and modelling errors，is indicated by the differences in the curves(De Jong van Lier and Wendroth，2016).Point 14 was visually identif ied as the closest to provide mean results for the cumulative functions，manifesting the average behavior of soil structure as a“point sensor”for monitoring SPQ in this f ield.

Fig.2 Soil water retention curves at the 15 sampling points for the assessed subtropical soil under maize crop cultivation in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil.The curve for location 14(bold solid line)was closest to the mean for all 15 curves.

TABLE II Results of the parameters a)of the Genuchten-Mualem model(see Eq.1 for this model)at the 15 sampling points for the assessed subtropical soil under maize crop cultivation in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil

Basic properties and SWRC-based indicators

As identif ied in many tropical and subtropical soils(Ottoniet al.，2019)，the soil texture was classif ied as clayey(Soil Survey Staff，2014)(Table I).Usually，clay intrinsic characteristics assign several soil properties，such as high microporosity(maximummicro=0.580 m3m-3)，large water retention，and low permeability(minimumKs=0.89 mm h-1);these features can be seen in Table I.However，the huge inf luence of the structure on these properties was evidenced at Points 9 and 15，where large values ofKsandmacrowere observed.The large SD forKsmeasurements at Points 9 and 15 should be taken into account.These results corroborate those of Dexteret al.(2004)and Sandinet al.(2017)，who observed that texture，Ks，and soil structure are strongly related.

Reynoldset al.(2008)emphasized that the ideal SPQ threshold that considers both ideal root zone water uptake and drainage is obtained when 18≤Ks≤180 mm h-1.The ranges 0.36≤Ks≤18 mmh-1and 180＜Ks≤360 mm h-1are considered intermediate thresholds of SPQ.Following these limits，SPQ at all the evaluated points was classif ied as ideal or intermediate.TheKsrange of 0.89—138.1 mm h-1，with larger SD for points whereKswas greater，indicated the strong variability of this property in the studied f ield.The large variability ofKsin space and time and its proper description are the theme of numerous discussion，including the choice of the constant or falling head method，f ield or laboratory measurements，and reaching full saturation and sampling according to adequate representative elementary volume(REV)(Laurenet al.，1988;Ottoniet al.，2019).Each decision implies an additional source of variation，whereas the choice of the relatively small number of three subsamples carries the risk of not totally encompassing the f ield variation.Although the mean values were taken from nine values(three cores×three replicates)，the large SDs indicate that the mean values should be cautiously examined before further use(Table I).

The measured mean OC content of 39.4 g kg-1(Table I)is within the ideal threshold(30—50 g kg-1)for the growth of most plant species，as pointed out by Reynoldset al.(2009).Only those of Points 5，9，and 13 were outside of this range.Point 9 presented the smallest value of OC(22.2 g kg-1).At this point，with its specif ic conditions，the largeKsvalue could be associated with the relatively small water retention.However，at Point 15，OC was 40 g kg-1，indicating that this measurement may not be sensitive enough to reveal soil structural changes.

Usually，ρb(Table I)indicates soil aeration，penetration resistance，water storage capacity，and permeability.Based on the study of Reynoldset al.(2009)，who described the threshold values ofρbbetween 1.25 and 1.30 kg dm-3as upper limit of good SPQ，crop growth may be inhibited in seven locations because of insufficient soil aeration and water f lux.The largest values ofρbwere observed at Points 10(1.40 kg dm-3)and 11(1.42 kg dm-3)，which could be associated with excessive mechanical resistance to root elongation，according to the same authors.However，the WRaindex had rather different values at Points 10 and 11(WRa=18.6 and 27.2 hPa m3m-3，respectively)，emphasizing differences in SWRC shapes(Table II).Because Points 10 and 11 exhibited similar texture andρb，WRaprovided additional information showing sensitivity for SPQ analysis(Armindo and Wendroth，2019).As already pointed out by McQueen and Shepherd(2002)，ρband penetration resistance may not be sufficient for an accurate assessment of crop relationship with soil structure;therefore，WRaappears to be an alternative to contribute to this analysis.

Following the threshold values for AC and PAW presented by Reynoldset al.(2009)，all points evaluated revealed optimum AC because their values were larger than 0.14 m3m-3.At the same time，nine of the 15 mean values of PAW(Table III)were smaller than 0.10 m3m-3and ergo classif ied as“poor”or“droughty”SPQ，regarding the ability of soil to store and provide water for plant roots.Castelliniet al.(2019a)and Cullottaet al.(2016)discussed the need of a superior limit for AC because it was not realistic to presume that the SPQ increases indef initely with soil aeration.They presented an ideal AC range between 0.10and 0.26 m3m-3for soils under agricultural and forestry use.According to this classif ication，only Points 4，5，and 6 presented ideal AC.The negative relationship between AC and PAW(Table III)indicates that when the soil has an AC of approximately 0.10 m3m-3，the value of PAWshould be approximately 0.3 m3m-3，indicating that a superior limit should be considered.

Point 13 presented the largest values of PAW(0.19 m3m-3)and WRa(135.3 hPa m3m-3)，which might be related to the high OCcontent(54.8 g kg-1).Points 2 and 3 displayed the same value for AC(0.29 m3m-3);however，the difference of 53%inAabetween these two points manifested their structural differentiation.The same was noticed at Points 4 and 7(PAW=0.14 m3m-3)，which had a difference of 67%in WRa.In both cases，different results ofθs，θ(1.0 m)，andθpwpgenerated the same values of AC and PAW，suggesting similar behavior in terms of SPQ.However，values ofAaand WRaref lected substantial differences in SPQ and soil structure，which were not detected by AC and PAW.Bacheret al.(2019)evaluated the performance of PAWas an SPQ indicator and found locations where PAWwas not sensitive enough to pick up differences between soil samples.However，these differences were captured by energy indices.Therefore，AWris considered a dynamic indicator capable of ref lecting temporal variations in SPQ with high precision.The lack of sensitivity of PAWand AC is related to the fact that they are parameters that def ine“total capacity”in terms ofθrange but do not reveal the arrangement of the soil matrix and pore connectivity and distribution.

TABLE III Soil physical quality indicators a)derived from soil water retention curve and pore size distribution functions at the 15 sampling points for the assessed subtropical soil under maize crop cultivation in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil

TABLE IV Pearson correlation coefficient(r)matrix for the measured soil properties and energy indices a)at the 15 sampling points for the assessed subtropical soil under maize crop cultivation in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil

The relationship between the mean measured soil properties and energy indices was evaluated based on the results of the Pearson correlation analysis(Table IV).No signif icant linear correlation between the relative hydraulic energy indices and any other physical properties was observed，evidencing the indices’contribution to SPQ analysis because their information was not captured by any other important and widely used key variables.Additionally，macropresented moderate correlations with the mean ofKs(r=0.634，P＜0.05)and withmicro(r=-0.592，P＜0.05).A negative correlation betweenKsandρb(Table IV)was also observed by Perkinset al.(2007)，who analyzed soils affected by mechanized military training and observed an increase inρbassociated with a decrease inKs.The energy indices presented very strong correlations when they were compared among themselves，with exception of that betweenAaand WRa，which was weak and not signif icant(r=-0.213，P=0.445).The absolute energy indices presented stronger correlations with measured physical properties than the relative ones;however，without signif icance at 5%in some cases.The indexAahad the strongest correlation with the measured properties，presenting moderate correlations with silt content(r=0.569，P＜0.05)，geometric mean ofKs(r=-0.462，P=0.083)，micro(r=0.608，P＜0.05)，andmacro(r=-0.517，P＜0.05).This result indicates the sensitivity ofAato capture the structural behavior of porous media because it is related to properties that are commonly used to refer to it.The index WRapresented a signif icant moderate correlation withρb(r=-0.541，P＜0.05)andmicro(r=0.532，P＜0.05).

The amount of energy associated with the gravitational potential to drain soil moisture fromθstoθfcis shown in Fig.3a.This hydraulic energy function reveals soil physicalstructural quality related to drainage after a rainfall event with low intensity and long duration.The soil hydraulic energy to store water over the PAW range for maize water uptake to promote its physiological functions is presented in Fig.3a.Because WRacorresponds to the total demanded energy betweenθfcandθpwp，when 50%of available water is consumed，i.e.，half of thexaxis in Fig.3a is reached，the energy used by the soil under a maize crop ranges from 11%to 35%of WRa.This result indicates the low energy by this soil to hold water in its pores(less than 5 hPa m3m-3related to 50%of normalized PAW).Armindo and Wendroth(2016)applied the energy indices and calculated AC and PAW for a tropical clay loam from Boca da Mata City.They found large AC and near zero PAWvalues and similar cumulative functions，indicating the possible existence of a pattern for soils from tropical and subtropical areas.This behavior was not observed by Bacheret al.(2019)，who applied the energy indices to Haplic Cambisols and Haplic Stagnosols collected in Ireland and found much larger values for WRathan those observed in the tropical soils.Huf dos Reiset al.(2019)applied theAaindex and the hydraulic energy function of aeration to evaluate soil structure of a Haplohumox under an integrated crop livestock system.Similar results for the cumulative energy function in the Cambisol evaluated in this work were observed for minimum tillage with winter crop use in the top layer.For conventional tillage，no-tillage，and chiseled no-tillage systems，the results forAa，and consequently，for the cumulative function were approximately 90%lower than those of minimum tillage.

Fig.3 Energy used by gravity to drain the normalized moisture in soil pore space from saturation to the f ield capacity point(a)and by soil to hold the normalized moisture in soil pore space from f ield capacity to the permanent wilting point(b)at the 15 sampling points for the assessed subtropical soil under maize crop cultivation in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil.The bold solid line(Point 14)best represented the mean curve.θs:saturated soil water content;θfc:f ield capacity;θpwp:permanent wilting point;θj:any soil water content betweenθs andθfc;θi:any soil water content betweenθfc andθpwp.

PSD-based indices

Based on the suggestion of Blott and Pye(2001)，who analyzed the particle size distribution curve，the assessed points in this work were classif ied with respect to their PSD curves.Most of the points f itted to“coarse-skewed”(skewness between-0.1 and-0.3)and“leptokurtic”(kurtosis between 1.11 and 1.5)classes.Kurtosis represents the degree of concentration relative to the mean value，indicating that the pore distribution revealed a more peaked distribution at the center and more tailed at the extremes than the lognormal curve.The case SD=1 indicates no variation in pore diameter and that pore diameters are approximately the same size.The values of SD ranged from 1.81 to 2.61，being smaller than those observed by Armindo and Wendroth(2016)for both tropical and temperate soils.Skewness represents the symmetry or preferential dispersion around the mean.The negative skewness values shown in Table III indicate an excess of small pores(relative to the lognormal distribution)，which is expected in clay soils(Blott and Pye，2001;Reynoldset al.，2009).Point 13 presented a combination of the largest variation in the range of pores and the greatest excess of small pores.

Several authors，such as Ghibertoet al.(2015)and Shahabet al.(2013)，have used the indices based on the PSD function proposed by Reynoldset al.(2009)to assess the structural quality of some soils.Comparable analyses were performed to indicate optimum ranges for clay loam soils(Typic Argiaquoll)collected in Canada(Reynoldset al.，2009)，silt loam and silt soils(Molisols)collected in Argentina(Ghibertoet al.，2015)，clay soils under different soil tillage systems collected in Italy(Castelliniet al.，2013)，and sandy to clay loam soils(Haplocalcids)collected in Iran(Shahabet al.，2013).The ideal range in each study was established based on a particular soil type，among all others，that had properties considered desirable regarding SPQ(Table V).The results of the mean indices for this subtropical soil can be found in Tables III.Comparing our results with the values found in the literature，only the mean value of kurtosis(1.13)(Table III)was within the range proposed by Reynoldset al.(2009)(Table V).The indicesdmode，dmedian，dmean，SD，and skewness were not within any proposed optimum range.Castelliniet al.(2014)evaluated soils cultivated with winter wheat and observed values outside the optimal values or reference curves proposed by Reynoldsetal.(2009).Nevertheless，for a tropical soil under native forest in Maringá，Brazil，Armindo and Wendroth(2016)founddmode，dmedian，anddmeansimilar to those presented in Table V，following the pattern observed for the energy indices.These results corroborate those of Cullottaet al.(2016)，Pulido-Moncadaet al.(2015)，Iovinoet al.(2016)，and Castelliniet al.(2013)，indicating that the optimum values proposed may not be appropriate to a wide range of agricultural soils;therefore，optimum ranges should be no more than indications regarding possible values.

Spatial variability of the soil properties and energy indices

The parameters of the semivariogram models that best f itted the experimental semivariograms are presented in Table VI.The best f itting was found using the spherical and Gaussian model，with practical ranges(a)between 14.3 and 131.8 m.As performed by da Silvaet al.(2017a)and Castelliniet al.(2019b)，who applied the criterion suggested by Cambardellaet al.(1994)，the DSD was evaluated based on the nugget-to-sill ratio.A strong spatial dependance was found for all semivariograms of hydraulic energy indices because they presented DSD＜25%，which is basically a result of the small nugget(C0)effect.Strong spatial structural behavior manifested by a low DSD was also observed for clay content，ρb，Ks，andmicro.The spatial dependance was moderate(DSD between 25%and 75%)for sand and silt contents and weak for OC(DSD＞75%).

TABLE V Desirable ranges of mode，median，and mean of pore diameter(d mode，d median，and d mean，respectively)，standard deviation(SD)，skewness，and kurtosis derived from pore size distribution function proposed in the literature

TABLE VI Parameters a)of the semivariogram models for the soil physical properties and for the hydraulic energy indices of soil physical quality

The maps of spatial variability of soil properties and indicesAa，WRa，Ar，WRr，and AWrare shown in Fig.4，presenting the spatial changes of SPQ indices throughout the f ield and ref lecting the spatial dependence of the structural information.An inverse relationship was obvious between the maps ofArand WRr，in which the regions with large values of WRrcorresponded to zones with small values ofAr，andvice versa.The map of AWrpresented the same spatial distribution pattern.This result was expected because all these three indices are related to effective pore space(θs-θpwp).Moreover，Arand WRrare normalized by the same energy proportion(Aa+WRa)and both had very strong correlations with AWr(Table IV).Interestingly，the spatial distribution pattern ofAawas different from those of the other indices，which can be explained by the higher values ofC1and VER.The maps of texture and other soil properties exhibited different spatial behavior，although some similarities can be observed between those ofρbandmicro，with large values ofρbassociated with smallmicroin the center.This relation was conf irmed by the results presented in Table IV，where a negative linear regression between these two attributes was found.

Fig.4 Spatial variability of soil physical properties and indices for the studied f ield in the Applied Agricultural Sciences Farm of the Federal University of Paraná，near Pinhais City，Paraná，Brazil.ρb:bulk density;K s:saturated hydraulic conductivity;OC:organic carbon;micro:microporosity;A a:absolute aeration energy;WR a:absolute water retention energy;A r:relative aeration energy;WR r:relative water retention energy;AW r:relative air-water retention.Geometric means were used for K s.The units for the properties and indices are as follows:%for sand，silt，clay，A r，WR r，and AW r;kg dm-3 forρb;mm h-1 for K s;g kg-1 for OC;m3 m-3 for micro;and hPa m3 m-3 for A a and WR a.

Because the range shows the length of spatial dependence，the smallest range forKs(a=14.3 m)suggests that this variable was spatially correlated only over short distances.For the energy indices，the range was also small and close to or smaller than the sampling grid(Table VI)，indicating that these properties could be better represented by a f iner grid.

Spatial dependence and small variation in the area was expected because the property analysis was conducted in a spatially continuous area，without sudden spatial changes.However，for properties with intrinsic suddenness，such asKs，the spatial variation seemed to not be fully represented，indicating that both REV and the sampling grid should be revised.

CONCLUSIONS

In this work，SPQ indicators based on measured and calculated soil physical properties were applied to examine a clayey subtropical soil under a tilled maize cropping system.The experimental results ofρb，macro，micro，φ，andKsexhibited similarity with soil properties found in other soils under agricultural use.The propertyKspresented large variability amongst subsamples and over the entire f ield and showed a signif icant correlation with OC content.The indicators based on PSD functions indicated an excess of micropores，whereasdmode，dmedian，dmean，and SDindices did not f it the optimum ranges established for other reported soils in the literature，indicating that those threshold values should be carefully examined before they are applied for a particular f ield site.Future studies might also focus on the effects of sample volume and REV on the resulting structural indices.The energy indices were useful in the analysis of sample points where commonly used properties and other indices presented similar behavior.Furthermore，hydraulic energy indices and functions indicated suitable structural conditions for soil aeration and water retention processes，and possible threshold values were described for the evaluated area.

ACKNOWLEDGEMENT

This study was f inanced in part by the Coordination of Improvement of Higher Education Personnel(CAPES)，Brazil(No.001).