A novel dowscaling procedure for compositional data in the Aitchison geometry with application to soil texture data

In this work, we present a novel downscaling procedure for compositional quantities based on the Aitchison geometry. The method is able to naturally consider compositional constraints, i.e. unit-sum and positivity. We show that the method can be used in a block sequential Gaussian simulation framework in order to assess the variability of downscaled quantities. Finally, to validate the method, we test it first in an idealized scenario and then apply it for the downscaling of digital soil maps on a more realistic case study. The digital soil maps for the realistic case study are obtained from SoilGrids, a system for automated soil mapping based on state-of-the-art spatial predictions methods.


Introduction
Uncertainty Quanfitifcation (UQ) is a crucial aspect for numerical tools intended to simulate physical processes, since it is important to provide an extensive analysis of the uncertainty of the outputs related to the variability of the inputs. Classical methods to perform this task are based on Monte Carlo (MC) simulations [22]. Here, an ensemble of realizations of the input parameters is used to feed a mathematical/numerical model, aiming to assess the distribution of the response in the face of uncertain inputs. In this broad framework, whenever parameters are characterized by a spatial distribution, geostatistical stochastic simulation can be employed to generate input scenarios for the model [6]. The geostatistical approach allows one to account for the spatial dependence characterizing the input parameters and to model the spatial structure expected for the realizations (range of variability, degree of smoothness) through a spatial covariance function. Nonetheless, sound geostatistical simulation needs to take into account the possible constraints of the data, particularly when these represent compositional information. For instance, soil moisture retention plays an important role in models that simulate hydrogeological processes and depends on a number of terrain properties, such as the soil texture. The latter in turn is determined by particle-size fractions (psfs), i.e. the relative percentages, in terms of soil composition, of clay, silt and sand, the three categories in which grains of fine earth are divided depending on their size, see e.g. [26]. When some sparse samples are available, geostatistical techniques such as Kriging and conditional Gaussian (co)simulation can be used to interpolate the available observations and assess the associated uncertainty. However, neglecting the inherent characteristics of these data may result in inappropriate results, such as prediction of negative components or modeling spurious correlations [23]. These serious limitations hinder the use of classical geostatistical methods based on the Euclidean geometry in the presence of compositional data (see, e.g. [1,7]).
In the last decades, an increasing attention has been devoted to developing analytics tools able to account for the features of compositional data, starting from the work of Aitchison [1]. Nowadays, Compositional Data Analysis (CoDa, [12,36]) is a wellestablished area of statistics, which studies models and methods for compositional data, grounded on the Aitchison geometry for the simplex. The Aitchison geometry is based on the foundational idea that, in compositional vectors, only the log-ratios among components represent a meaningful information to be accounted for in the statistical analysis. In a geostatistical setting, this foundational idea led to the development of new kriging methods for compositional data, which were successfully applied in several environmental studies [37,43].
In this work, we focus on the problem of geostatistical downscaling of compositional quantities. This is relevant in applications where no (or limited) direct observation is available within the study area -because of cost or environmental constraints -but low-resolution information is available across the region. This is the case of our motivating study, which focuses on the stochastic characterization of soil texture within a mountain river catchment, aiming to model the hydrogeological instability -and consequent natural hazard -of the region. In this case, no direct observation of particle-size fractions is available, but low-resolution data are reported in public databases, such as SoilGrids [19,21]. In this case, characterizing the spatial distribution of the soil texture requires to operate a change of support of the available (compositional) information and to assess the corresponding uncertainty.
To the authors' knowledge, none of the available methods for (geostatistical) downscaling allows to account for compositional constraints. For instance, methods of areato-point kriging and stochastic simulation available in the literature [24] inevitably suffer from the limitations of the Euclidean methods. We here propose an extention of Area-To-Point Regression CoKriging (ATPRCoK) -and associated stochastic simulation -to compositional vectors that, based on the Aitchison geometry, allows to overcome such issues and provide stochastic scenarios for the target compositional parameters.
The remaining part of this work is organized as follows. In Section 2 we recall the area-to-point regression (co)kriging method; in Section 3, we present the downscaling prediction framework for compositional data; in Section 4 we recall the definition of psfs, which are used in Section 5 to exemplify and test the features of the method in a first to synthetic case and then in a real scenario. Finally, we apply the method to a case study within the Caldone catchment in the Northern Italy city of Lecco, where we show how the method is able to provide psfs data at a length-scale most of the time very difficult or impossible to be determined.

Area-to-point regression kriging
In this section, we recall the main features of Area-To-Point Regression Kriging (AT-PRK) and Area-To-Point Regression Cokriging (ATPRCoK); for further details see e.g. [46,47]. Let us consider a scalar continuous random field {Z(x), x ∈ D} defined over a geographical region D ⊂ R d . Let us discretize Z(x) as where Z j denotes the discretized element at pixel j, ν j defines the geographical support of the j-th pixel having center x j ∈ D, and |ν j | denotes the measure of the support ν j . We assume the measure of the pixel support to be equal for all the pixels covering the region D and consider two levels of spatial resolution, one coming from a coarse discretization, denoted by the index K = 1, . . . , M , and another coming from a fine discretization, denoted by the index k = 1, . . . , N . The fine and coarse supports are such that we can define an integer number P = |ν K | |ν k | . Moreover, when using a Euclidean geometry for the data -which is the standard setting for which ATRPK is developedthe low-resolution random field is assumed to be obtained as an arithmetic mean of the high-resolution one, i.e., for K = 1, . . . , M, Starting from one complete realization of the low-resolution field Z K , we want to estimate the high-resolution field Z k , i.e. perform downscaling. ATPRK allows one to compute an estimate of the field Z k as a (linear) combination of two parts: regression and Area-To-Point Kriging (ATPK), see e.g. [3,15,24,25]. It uses a linear regression model on a set of covariates for the mean term of Z k , and kriging to interpolate the residuals from the regression model. The ATPRK predictor Z k of the field Z k at a given fine scale pixel ν k is defined as with β l , l = 1, . . . , L and λ K , K = 1, . . . , M unknown real quantities. The first sum in (3) is a classical linear regression term, describing the mean of Z k where u l k are a set of known fine-resolution regressors. Given a realization of the field Z K , one may linearly upscale Equation (4) to obtain a linear regression model for Z K , i.e., where u l K are the upscaled regressors, i.e.
By combining equations (2)-(5)-(6), the regression coefficents β l can be thus estimated by using a standard fitting procedure (e.g. Ordinary Least Squares (OLS) method) on the low-resolution field Z K , see e.g. [20,28]. The second term in Equation (3) is the Area-To-Point-Kriging (ATPK) term. It is the best linear unbiased predictor from the coarse residuals e K , defined as e K = Z K − E[Z K ]. In ATPK, the residual at a given fine pixel k is predicted as the best linear combination of the coarse residuals, subject to unbiasedness, i.e., e k = K λ K e K , where e k is the fine resolution predicted residual and λ K solve In practice, the ATPK predictor is often computed from a subset of M < M of residuals (typically selected in a neighborhood of the target pixel), to reduce the computational burden of the procedure. The optimal weights λ = (λ 1 , . . . , λ M ) are computed by minimizing the prediction error variance, which yields the following kriging linear system Σ 1 4 Here, the element in position (K 1 , K 2 ) of matrix Σ is the block-block covariance between the residuals at the coarse pixels centered at x K 1 and x K 2 ; the K-th element of σ is the point-block covariance of the residuals between fine and coarse pixels respectively centered at x k and x K , and µ is a Lagrange multiplier. Note that, in practice, neither the residuals nor their covariance are observed, but need to be estimated from the data. Residuals are typically estimated by difference from the estimated regression term. Estimating the covariance structure is more critical. Under the assumption that the residual field e k is stationary and isotropic and denoting with C k 1 ,k 2 the covariance between pixels k 1 and k 2 of the residual at fine scale, we can compute the block-block covariance as The point-block covariance is then given by The critical point of the ATPK method is thus the determination of the covariance structure at the fine scale, which cannot be directly estimated, as the data are given at the coarse scale only. The problem of estimating the fine-scale semi-variogram γ k 1 ,k 2 from coarse-scale data is known as a deconvolution problem. We shall not focus on this problem in this work, as it is completely analogous to that encountered in the Euclidean setting -for which a number of methods are available. We refer to [15] for more details on the deconvolution method used in this work. In case of multi-dimensional random fields, the ATPRK framework changes slightly in order to take into account possible cross-correlations among field components. Generalization of ATPRK to the multivariate case is analogous to cokriging, and yields Area-To-Point-Regression-CoKriging (ATPRCoK), see e.g. [47]. In ATPRCoK the coarse residuals appearing in (3) are replaced by the residuals of all the components of the multi-dimensional random field, in order to consider possible cross-correlations among their components. If we consider a p-dimensional random field {Z(x), x ∈ D}, its ATPRCoK discrete prediction is, where β l ∈ R p are the vectors of the unknown regression coefficients. The matrix Λ K ∈ R p×p contains the unknown cokriging coefficients and e K ∈ R p is the vector of the coarse scale residuals. The optimal weights Λ K are found by solving a system analogous to (8), but considering covariances and cross-covariances within/among fields components, as in a standard cokriging setting.

5
In this section, we consider the problem of downscaling compositional data and we propose a method which extends ATPRCoK to the Aitchison geometry, and naturally takes into account the compositional nature of the data.

Compositional data in the Aitchison simplex
A compositional data point Z = (Z 1 , ..., Z p ) is typically represented as a vector whose elements are proportions (or percentages) of a whole, named total. In this case, compositional vectors are characterized by the unit-sum constraint i Z i = 1, where we denote with Z i ≥ 0 the i-th component of compositional data point. More generally, compositional vectors are data which convey only relative information, being subject to a constant-sum constraint. Here the total is typically of no interest for the analysis, in the sense that expressing the data w.r.t. a different total (i.e., in proportion, percentages or ppm) should not change the results of the analysis (i.e., scale invariance).
Because of the range limitation and the possible spurious correlation of compositional vectors [23,30], the Euclidean-based statistical framework was proved to be ineffective for the spatial prediction of this type of data, although a number of authors have ignored this aspect, see e.g. [10]. Other works (e.g. [45]) tried to account for the particular nature of regionalised variables expressing relative fractions by proposing an extension of kriging called Compositional Kriging (CK). CK predictions respect the constraints of positivity and constant sum value. However, the CK algorithm is based on empirical considerations rather than a coherent probabilistic model, and is therefore not suited for stochastic simulation. Our developments follow the direction of research on compositional kriging explored by [31,44], who formulated geostatistical models and methods based on the Aitchison geometry for the simplex (see [35,43] for recent reviews). Presently, the standard approach to the statistical analysis of compositional data is the one pioneered by Aitchison, see [1], which is based on the particular geometry of the simplex [2,36,33]. A p-dimensional compositional vector Z = (Z 1 , . . . , Z p ) is an element of the p-dimensional standard simplex, S p , which is defined as In [2,33] group operations are defined to give the simplex a structure of a real vector space. These are the perturbation ⊕ (sum) and powering (product by a constant) operations, defined, for x, y ∈ S p , and α ∈ R, respectively as Here, C(·) denotes the closure operation 6 The space S p can be equipped with a (finite-dimensional) Hilbert space structure when considering the Aitchison inner product, defined, for x, y ∈ S p as The defined inner product induces a norm · a := ·, · a , which in turn induces a distance d a (x, y) = x y a , x, y ∈ S p , where x y denotes the perturbation of x with the reciprocal of y, i.e., x y = x ⊕ ((−1) y). The Hilbert space structure identified by these operations is called Aitchison geometry, or Aitchison simplex [33].

ATPRCoK in the Aitchison geometry
The statistical approach proposed by Aitchison [1] and following authors [36,34] consists in analyzing compositional data in the context of the Aitchison geometry. Here, a large number of multivariate statistical methods grounded on a geometric perspective (e.g., principal component analysis. regression) can be properly reformulated to account for the inherent properties of compositional data. From an operation viewpoint, the standard procedure of analysis consists in transforming the original data by applying an isomorphism from the p-dimensional Aitchison simplex to the classical Euclidean space R p−1 or R p , perform the statistical analysis on the transformed data and finally back-transform the results in the original space. This strategy was proved to be fully equivalent to working directly in the Aitchison simplex for a number of statistical methods (see, e.g. [12]). In this section, we shall formulate ATPRCoK in the Aitchison simplex, and then prove that one may equivalently perform the analysis by relying on the so-called Isometric Log-Ratio (ILR) transformation, which is an isometry that maps the simplex in R p−1 . The latter associates to a compositional vector z ∈ S p the coordinates of this vector with respect to an orthonormal basis of the simplex (ψ 1 , . . . , ψ p−1 ), i.e., Note that the ILR is a linear transformation, see e.g. [12]. For several compositional methods (e.g., principal component analysis, regression, see, e.g., [36]), it was shown that the choice of the basis does not influence the results of the analysis. However, specific choices for the basis can lead to practical advantages. For instance, the basis could be chosen in such a way as to grant uncorrelation of the resulting transformed data, or to ease the interpretation of the results (see, e.g. [13]). In [11] the authors used ATPCoK to downscale psfs data transformed with Additive-Log Ratio (ALR), using the silt fraction as a reference for the ratios. ALR was there used as a practical solution to account for the compositional nature of the data, but the modelling assumptions were not explicitly stated, and the results were interpreted only in terms of prediction accuracy with respect to a given test set. Recent works highlighted some limitation of the ALR transformation, as this is not isometric, thus does not preserve the Aitchison geometry (see, e.g., [36]). In this section, we propose a general method for the statistical downscaling and simulation of compositional data which extends the ATPRCoK to the context of the Aitchison simplex. We call the method ILR-ATPRCoK to recall the computational strategy we propose to perform ATPRCoK in the Aitchison simplex, which is based on the ILR transformation. Here, we shall also prove that the strategy based on ILR is fully equivalent to working directly in the simplex itself.
In the following, we denote by {Z(x), x ∈ D} a random field valued in S p , defined over a Euclidean region D ⊂ R d . To indicate the Aitchison center of the field (i.e. the mean value in Aitchison geometry), the Aitchison covariance and the integral operator of Z(x) over a region τ ⊂ R d , we use respectively the same notation used in [27,32], that is • Aitchison covariance operator, acting on a (non-random) element z ∈ S p , for We refer to [2,12] for an insight of the geometry of the random compositions in the Aitchison simplex, and to [5] for a recall on covariance operators in Hilbert spaces. For the element Z(x) of the compositional field at x ∈ D, we assume the following model with e(x) the residual term. We model the center of the field through a linear model in where β l ∈ S p , l = 1, . . . , L are the vectors of unknown regression coefficients and u l (x) ∈ R, x ∈ D, are the covariates. Furthermore, we assume that the residual is second-order stationary, see [43], i.e., that the covariance structure (in Aitchison geometry) for a random composition Z(x), x ∈ D, only depends on the increment among locations To simplify the notation, we shall indicate the stationary covariance function C a simply by C a .

8
With a notation analogue to that used in Section 2 we consider the discretized versions of the field {Z(x), x ∈ D}, denoted by Z k (resp. Z K ) and obtained at a fine (resp. coarse) discretization scale, namely Here, the powering by 1 |ν k | and 1 |ν K | is intended as acting element-wise. Given the realization of the coarse-scale field Z K , and by analogy with (3), we define the ILR-ATPRCoK predictor of the fine-scale field Z k ∈ S p as where Λ K ∈ R p×p is a matrix of ATPCoK unknown weights to be optimized, is the matrix-by-composition multiplication -consistent with perturbation and powering, as defined in [36] (p. 55), i.e., denoting by e K,i , i = 1, . . . , p the elements of e K ∈ S p and by λ K,i,j , i, j = 1, . . . , p the elements of Λ K The residuals e K ∈ S p represent the upscaled residuals of (21), defined as Under the assumption that the regression coefficients β l and the covariance function C a are known, the optimal weights Λ k in (25) are found as to guarantee that the ILR-ATPRCoK is the Best Linear Unbiased Predictor (BLUP) in S p , i.e., by solving the constrained minimization problem arg min where µ is the spatially constant residual mean.
The following result states that, by applying the ILR transformation (i.e. mapping from the Aitchison simplex to an Euclidean space, via an isometric isomorphism), one obtains an equivalent formulation of problem (27) that can be solved using the standard ATPRCoK presented in Section 2. The proof of Proposition 1 is given in Appendix A.
Proposition 1 Given a compositional random field Z(x) valued in S p and a random field Y(x) valued in R p−1 defined as Y(x) = ILR(Z(x)) for x ∈ D, the BLUP in S p for Z k -found by solving (27) -coincides with the ILR-back-transformed ATPRCoK predictor for Y k defined in (3), i.e., Even though β l is rarely a priori known, an estimate of β l can be obtained by backtransforming the corresponding estimate of the coefficient vectors β l Y referred to the ILR-transformed field Y(x) (see, e.g., [36]). Similarly, an estimate of the covariance operator C a can be obtained from the estimated (Euclidean) covariance operator C Y of the vector field Y(x). In this work, for β l Y we shall consider OLS estimates, whereas for C Y the estimates obtained by Goovaerts' deconvolution [15] of classical crossvariograms.
Note that the equivalence between the ATPRCoK in the Aitchison simplex and the Euclidean ATPRCoK on ILR-transformed data (as stated in Prop. 1), implies the possibility to analogously perform Block Sequential Gaussian Simulation (BSGS) [42], as BSGS grounds on the same hypothesis as ATPRCoK, and it is indeed based on the latter method. In the context of Uncertainty Quantification (UQ), BSGS is key to propagate the uncertainty in numerical models that take as input downscaled compositional data, as we discuss in Section 6.
Finally, one should note that, since the assumptions are made with respect to the Aitchison geometry, the mass-preserving property as stated in the Euclidean framework, i.e. (see [24]), does not hold. In a discrete prediction setting, as in Section 2, the Aitchison geometry predictions respect the following centre-preserving property Indeed, since P = |ν K | |ν k | ∈ Z + , as defined in Section 2, one has This means that, in the Aitchison simplex, coarse areal data coincide with the geometric mean of the predicted fine areal values (normalized to having unit-sum).
In the following sections we exemplify the proposed methodology through its application to particle-size fractions, whose definition is recalled in the next Section 4.

Particle size fractions
Soil texture is a classification instrument used to determine soil classes. More specifically, soil texture is quantitatively determined on the basis of the relative fractions of the fine particles of different sizes that compose the terrain. Soil particles under 2 mm are divided in three groups • clay: particles with a diameter less than 2 µm; • silt: particles with a diameter comprised between 2 µm and 50 µm; • sand: particles with a diameter comprised between 50 µm and 2 mm.
Fractions of clay, silt and sand are usually indicated as particle-size fractions (psfs). Soil texture classes are determined by the relative percentages of psfs, according to a standard that may vary depending on the country. The most common classification is Figure 1: Soil texture triangle: Soil texture classification according to the USDA classification system, based on relative fractions of clay, silt and sand. Figure taken from [17]. that used by the United States Department of Agriculture (USDA) [41], which distinguishes twelve major soil texture classes shown in Figure 1. The classes are typically named after the primary constituent particle-size or a combination of the most abundant particles sizes, e.g. sandy clay or silty clay. A fourth term, "loam", is used to describe equal proportions of sand, silt, and clay in a soil sample, and leads to the naming of even more classes, e.g. clay loam or silt loam.

Validation
The geostatistical method outlined in the previous sections has been implemented in R-3.6 [39] using the libraries gstat [18,38] for ATPK and geostatistical simulation and compositions [4] for the analysis of compositional data. In particular, for the variogram deconvolution we use the Goovaerts' procedure [15,16]. We define a continuous ran- In our analysis, compositional data are transformed using the function ILR of the package compositions, see e.g. [4]. The basis used for the transformation is the one introduced in the original article [12], based on the partition of the vector of compositional variables in two sub-compositions, the first consisting of Z 1 and Z 2 and the second containing only Z 3 . For the sake of illustration and in view of the motivating study, we shall interpret Z(x) as the psfs at x (i.e., Z 1 (x), Z 2 (x), Z 3 (x) represent the composition in clay, silt and sand, respectively). Nonetheless, the validity of the simulation study here presented is clearly not limited to the specific case of psfs.

Synthetic data
To assess the performance of the proposed method, we consider a simulated dataset Z(x), x ∈ D, with support measure |ν k | = 20 × 20 m 2 , on a given rectangular domain D with area |D| = 10000 × 9160 m 2 . The compositional vector Z(x) is modelled as a process with a given spatially-constant center Cen(Z(x)) = µ and stationary-isotropic covariance structure. From the operational viewpoint, the mean µ = C[(µ 1 , µ 2 , µ 3 ) ] is set based on independent uniform distributions µ i ∼ U [0, 1], i = 1, 2, 3. Compositions were simulated by back-transforming through ILR −1 two-dimensional Gaussian random vectors Y, with constant mean µ Y = ILR(µ), and stationary-isotropic marginal variograms from the spherical model without nugget [8,9]. For the following simulations, the components of Y are always assumed to be uncorrelated, and the marginal ranges are both set to 2000 m. In each simulation, the common sill σ 2 is sampled according to a uniform distribution U [0.025, 2.5]. In Figure 2 we show an example of realization of the psfs distribution.
Starting from this set of synthetic psfs, we perform a sequence of upscaling-downscaling procedures, as follows. Downscaling is done using either ATPRCoK or ILR-ATPRCoK and upscaling either in Euclidean or Aitchison geometry, so that four different possibilities arise. In the following, we call AA the upscaling in the Aitchison simplex and downscaling via ILR-ATPRCoK, EE the upscaling in the Euclidean space and downscaling via ATPRCoK whereas EA, AE are the mixed methods, referring respectively to upscaling in the Euclidean space and to downscaling via ILR-ATPRCoK and upscaling in the Aitchison geometry and downscaling via ATPRCoK.
For each method, we consider a set of 100 realizations of the fine scale compositional field, each yielding a reconstructed field after the upscaling-downscaling process. The upscaling factor P is set each time by randomly and independently sampling in the discrete range {2 2 , 3 2 , ..., 30 2 }. For each method and each realization we compute the   Figure 3: Synthetic data results: a) boxplots of the sample mean error. Blue points and red segments are respectively the spatial mean and spatial standard deviation of the sample mean error. b) the mean (across realizations) number of pixels that violate the positivity and unit-sum constraint for the four methods considered. c) relative violation of the unit-sum constraint. A value of 1% on the x-axis indicates that the reconstructed psfs sums e.g. to 1.01. 14 sample mean error, i.e. the average of the Euclidean distance between initial and reconstructed psfs fields -the average being taken over the realizations. The sample mean error is not computed through the Aitchison distance as this is not defined for recontructed psfs violating the compositional constraints.
Even if the distribution of the sample mean error between the considered methods, reported in Figure 3 a), would suggest a substantial equivalence among the methods ,  Figures 3 b), c) clearly show that, unlike ATPRCoK, ILR-ATPRCoK is able to produce psfs maps that are consistent with the unit-sum and positivity constraints. Indeed, the ATPRCoK shows a violation of the positivity constraint in about 10 3 pixels on average, representing roughly 0.5% of the study area. We then perform a sensitivity analysis with respect to random perturbations of the initial data of the downscaling procedure, for the four methods described above. This case is representative of input data characterized by a given degree of uncertainty. We thus consider a realization of the synthetic psfs field in case of µ = ( 1 3 , 1 3 , 1 3 ) and sill σ 2 = 0.1 (i.e. as in Figure 2) and we set the upscaling factor to P = 225. Let us indicate with K = 1, . . . , M the elements of the coarse maps. The upscaled maps Z K are then perturbed with a set of i.i.d. Gaussian random errors K . Similarly as before, these perturbations were generated on the ILR transforms, by adding a zeromean independent Gaussian error with variance s 2 , ranging from 10% to 100% of the sill σ 2 .
In Figure 4 a), b) we report the boxplots of the error maps for each value of s 2 . The error maps are computed, for each pixel, as the Euclidean distance between initial and predicted psfs. We note that both ATPRCoK and ILR-ATPRCoK are quite robust even in case of relatively high perturbations of the initial data. In Figure 4 c), d) we show the histograms of the maximum, across simulations, of the violation maps of the unit-sum constraint. For instance, the vertical bar in correspondence of the range [1,1.2] in Figure 4 d) indicates the count of pixels whose maximum discrepancy (across simulations) from unity of the sum of psfs is between 1% and 1.2% (i.e., the sum is in [1.01,1.012] or [0.988,0.99]). These results clearly show the ability of ILR-ATPRCoK method to produce results consistent with the unit-sum constraint, unlike the ATPRCoK method which yields maps with a significant violation of the aforementioned constraint.

Downscaling SoilGrids data
In this section, we test the performances of the proposed method in downscaling psfs from soil digital maps publicly available. This case is considered to analyse compositional random fields having a realistic spatial distribution. For this purpose and in view of our case study, we consider SoilGrids, which is a system for automated digital soil mapping based on state-of-the-art spatial predictions methods [19,21] released in 2014 by ISRIC (International Soil Reference and Information Centre) -World Soil Information, a non-profit organization funded by the Dutch government. SoilGrids predictions are based on globally fitted models using soil profile and environmental covariate data. When first released, SoilGrids provided a collection of soil properties and class maps of the world at 1 km spatial resolutions, produced using automated soil mapping based on statistical regression models. In 2017, the resolution has been increased to 250 m and the accuracy of the predictions has been greatly improved by using machine learning algorithms instead of the previously employed linear regression [21]. In 2020, SoilGrids released a version where, among other updates, soil map predictions are provided with a mean value together with an uncertainty level map. SoilGrids data are available publicly under the Open DataBase License. Among SoilGrids predicted variables, relevant to this work are clay, silt and sand percentages at different soil depths. In this section, the values considered for geostatistical downscaling are those referred to the topsoil, i.e. depth of 0 cm.
We focus on a geographical domain with area |D| = 15750 × 16000 m 2 , located in a pre-Alpine area, more specifically the basin of Pioverna river in the Lombardy region in Northern Italy near the city of Lecco. This region was selected as it is similar, from the geomorphological viewpoint, to the area analyzed in the case study presented in Section 6. The psfs as available in SoilGrids are reported in Figure 5. Based on these data, we test the performance of the ILR-ATPRCoK method, at different levels of the upscaling factor P . Following the procedure described in Subsection 5.1, we consider a sequence of upscaling-downscaling operations, both in Aitchison and Euclidean geometry, of the SoilGrids data in Figure 5. For each upscaling factor P in the range {2 2 , 3 2 , ..., 10 2 } and for each pixel in D, we compute the Euclidean distance of the psfs estimates from the initial SoilGrids data, yielding a set of error maps (one for each value of P ). These are displayed through boxplots in Figure 6 a), b). We note that, mainly at high uspcaling factors, the ILR-ATPRCoK method shows slightly better behaviour with respect to the classical ATPRCoK, producing solutions that are closer to the reference ones w.r.t the results produced via ATPRCoK downscaling.
In Figures 6 c), d), we report the histograms of the maximum of the violation maps of the unit-sum constraint experienced during the set of upscaling factors, for the four cases being considered. Interpretation of these histograms is fully analogous to that in Figure 4. These results confirm that the ILR-ATPRCoK method is able to produce downscaled maps that are consistent with the unit-sum constraint, as opposed to the ATPRCoK downscaling method. Finally, we do not report any violation of the positivity constraint for ATPRCoK, differently from what shown in the tests reported in Section 5.1.

A case study
Our case study considers an application to a domain D centered on the city of Lecco, located in the Lombardy region in Northern Italy, which is crossed by three streams (Bione, Caldone, and Gerenzone) that have the typical characteristics of the pre-Alpine area. The hydrographic basin of the Caldone water course is 24 km 2 wide, with an altitude ranging from 197 m a.m.s.l. to 2118 m a.m.s.l. at the top of Grigna Meridionale mountain. Geologically, the basin is characterized by rocky outcrops in the higher part (mainly limestone and clastic rock), while downstream towards the city the river flows through a floodplain. The average precipitation over the city of Lecco is about 1400 mm/yr.
The combination between a short hydrologic response time, high slope, intense sediment transport and flow within a densely urban area makes the Caldone river a suitable case study for hydrogeological instability and hazard. This motivates the development of numerical models intended to simulate hydrogeological processes, such as the SMART-SED simulation tool (Sustainable Management of sediment transpoRT in responSE to climate change conDitions) [14] which is able to simulate sediment transport resulting from slope erosion. These models typically need to be initialized with psfs maps, with a resolution consistent with the Digital Terrain Model (DTM), to be able to model properly the hydrological processes taking place in the study region. However, field measurements of psfs are not available at the study site, which motivates  the use of public repository to obtain indirect information on these input data. SoilGrids psfs at the study region have a pixel support with measure |ν K | = 250 × 250 m 2 . In terms of the USDA classification, the soil texture of the SoilGrids data for the present case study falls into the loam category. This kind of soil texture, according to [40], is classified as fairly permeable soil with moderate infiltration rates and moderate runoff potential. In the following, these coarse-scale data are downscaled to the resolution of the DTM employed for the SMART-SED model, i.e. 5 m, using ILR-ATPRCoK, following the methodology described in Section 3. Together with ILR-ATPRCoK results, we here aim to provide random realizations of the psfs fieldsobtained via Block Sequential Gaussian Simulation (BSGS) -as demonstration of the ability of the method to produce stochastic compositional maps compatible with coarse scale data.
Based on SoilGrids psfs data, we define the following coarse resolution maps In the ILR-ATPRCoK model, we consider as covariates u l k , l = 1, . . . , L, the DTM and its square, driven by the parabolic relation displayed in the scatterplot in Figure 9. For the fine map predictions ILR 1,k , ILR 2,k we thus consider the model  The fitted values are plotted against the observed values in Figure 10.
To perform ILR-ATPK, the spatial correlation structure of the fine residuals e 1,k and e 2,k is estimated by applying the Goovaerts deconvolution procedure to the variograms fitted to the coarse residuals e 1,K and e 2,K (based on a spherical model with nugget), and by assuming e 1,k and e 2,k to be uncorrelated. The latter assumption is supported by the residuals' analysis (not reported for brevity). Once the fine variograms of the residuals have been estimated, it is possible to solve the ATPK linear system, according to (8). The downscaled ILR are then backtransformed in the Aitchison space in order to get downscaled psfs, see Figure 12, left column.
Finally, in Figure 12, right column, we show a sample realization for the downscaled psfs, obtained via BSGS. These stochastic maps shall form the cornerstone to evaluate the propagation of the uncertainty associated with the psfs through the SMART-SED model, and eventually assess the sensitivity of the sediment transport model to this information.

Conclusions
We have presented a novel downscaling method for compositional data, based on the ATPRCoK method in the Aitchison geometry, with application to the geostatistical downscaling of psfs data. We have tested the method first in the case of synthetic data and then on a dataset from the SoilGrids online repository. In particular, we have shown the ability of the method to automatically handle the compositional nature of the considered data. Indeed, the proposed method produces maps that respect the unit-sum and positivity constraints, as opposed to the classical ATPRCoK method that produces maps which are not consistent with the compositional constraints.
Validation on both synthetic and SoilGrids data show good performances of the method in downscaling, as well as robustness to the uncertainty of the input data. This is critical to the use of data from public repositories in local analyses, when point observations are not available, as they are naturally prone to uncertainty at a fine scale. While a full account of SoilGrids uncertainty will be the scope of future work, a relevant feature of ILR-ATPRCoK method -similarly as ATPRCoK in the Euclidean setting -is the possibility to easily incorporate point observations collected at the site, thus anchoring the downscaled maps (either kriged or simulated) to such observations [29]. For instance, at the time of writing, a campaign of data acquisition is under way in the Caldone basin, and will support the definition of (possibly improved) random psfs maps, to be used as input to the SMART-SED model discussed in Section 6.

A Proof of ILR-ATPRCoK proposition
In the following we propose a proof of Prop. 1, i.e. the equivalence of the ILR-ATPRCoK predictor (in the simplex S p ) to the classical ATPRCoK predictor (in the space R p−1 ) by applying an isometric isomorphism. The equivalence must be intended in the predictor, unbiasedness and optimality conditions. In the following we make extensive use of the ILR properties defined in [36] (p. [37][38][39][40][41][42][43] and the fact that ILR : S p → R p−1 extract the Fourier coordinates of a basis projection for the vector z ∈ S p , i.e., where the rows of Ψ = [ψ i,j ] j=1...,p i=1,...,p−1 are (compositional) vectors identifying an orthonormal basis of the simplex {ψ 1 , . . . , ψ p−1 } and y = [y i ] i=1,...,p−1 ∈ R p−1 is the vector of coordinates (i.e. of the Fourier coefficients) of the identified basis of the simplex.
Let us start with the predictor, applying the ILR to the ATPRCoK predictor defined in the Aitchison space S p (25), we get where β l = ILR −1 (β l Y ) and e K = ((e Y K ) Ψ) . Being e Y K,i , the i-th element of the vector e Y K , i = 1, ..., p − 1, we have In this way, we obtain the ATPRCoK predictor in the Euclidean space R p−1 , i.e.,