Soil Salinity Prediction and Its Severity Mapping Using a Suitable Interpolation Method on Data Collected by Electromagnetic Induction Method

: Salt mining and shrimp farming have been practiced in the Non Thai district and the surrounding areas for more than 30 years, creating saline soil problems. To solve the soil salinity problem, soil salinity prediction and mapping utilizing the electromagnetic induction method (EMI) and spatial interpolation methods were examined in the Non Thai district, Nakhon Ratchasima province, Thailand. The research objectives were (1) to predict soil salinity using spatial interpolation methods and (2) to identify a suitable spatial interpolation method for soil salinity severity mapping. The research methodology consisted of ﬁve steps: apparent electrical conductivity (ECa) measurement using an electromagnetic induction (EMI) method; in situ soil sample collection and electrical conductivity of the saturated soil paste extract (ECe) measurement; soil electrical conductivity estimation using linear regression analysis (LRA); soil salinity prediction and accuracy assessment; and soil salinity severity classiﬁcation and overlay analysis with relevant data. The result of LRA showed a strong positive relationship between ECe and ECa. The correlation coefﬁcient (R) values of a horizontal measuring mode (HH) and a vertical measuring mode (VV) were 0.873 to 0.861, respectively. Four selected interpolation methods—Inverse Distance Weighting (IDW), Ordinary Kriging (OK), Ordinary CoKriging (OCK) with soil moisture content, and Regression Kriging (RK) without covariable factor—provided slightly different patterns of soil salinity prediction with HH and VV modes. The mean values of the ECe prediction from the four methods at the district level varied from 2156.02 to 2293.25 mS/m for HH mode and from 2377.38 to 2401.41 mS/m for VV mode. Based on the accuracy assessment with the rank-sum technique, the OCK is a suitable interpolation method for soil salinity prediction for HH mode. At the same time, the IDW is suitable for soil salinity prediction for the VV mode. The dominant soil salinity severity classes of the two measuring modes using suitable spatial interpolation methods were strongly and very strongly saline. Consequently, the developed research methodology can be applied to conduct soil salinity surveys to reduce costs and save time in other areas by government agencies in Thailand. Nevertheless, to apply the EMI method for soil salinity survey, the users should understand the principle of EMI and how to calibrate and operate the EM device properly for accurate ECa measurement.


Introduction
According to the 2021 Food and Agriculture Organization (FAO) report, the total global area of salt-affected soils, including saline and sodic soils, is 833 million hectares [1]. These areas extend over Africa, Asia, Australasia, and America. The distribution of saline and sodic soils based on the FAO/UNESCO Soil Map of the World covers 397 million hectares and 434 million hectares, respectively. Most salt-affected land lies in arid and semiarid environments. Soil salinization is a significant threat to many countries' agricultural production systems and food security. Table 1. Distribution of saline soil and its severity in Nakhon Ratchasima province [14].

Distribution of Saline Soil and Its Severity Area (km 2 ) Percent
Strongly saline, salt surface ground more than 50% In this study, we introduce an electromagnetic induction (EMI) method, which has been applied and suggested by many researchers [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], to measure the apparent soil electrical conductivity for predicting and mapping soil salinity with a suitable interpolation method. The EMI method measures the electrical conductivity (EC) of undisturbed soil in the field. The reading responds to salt content, moisture, and clay levels. When calibrated and interpreted with soil test results, an electromagnetic (EM) sensor can provide helpful information on salt and water movement [31]. An EM device has two wire coils at each end of the cable. One coil transmits a low-frequency radio signal to the ground, while the second coil (the receiver) acquires the secondary (return) signal from the ground (Figure 1). Highly electromagnetically conductive ground (such as saline, clay, or wet ground) produces a stronger secondary signal than less conductive ground. By measuring the difference in strength between the original signal and the secondary signal, the EM device can measure the ground's conductivity [32]. At a national level, Khongnawang et al. [27] used an EM device, EM38, with field survey data and inversion software to map the distribution of clay and cation exchange capacity (CEC) in sandy and infertile soil in Ban Haet district, Khon Kaen, Thailand. In their study, they firstly examined a linear regression relationship between EM38 ECa with either topsoil (0-0.3 m), subsurface (0.3-0.6 m), or subsoil (0.6-0.9 m) clay and CEC. Then, the predicted equations were applied to map clay and CEC on three different levels.
Meanwhile, at an international level, Yao and Yang [34] used a mobile electromagnetic induction system, including Geonics EM38 and EM31, to perform an electromagnetic field survey and quantify the spatial pattern of soil salinity in Dongyihe Village of Kenli County, Shandong Province, China. In their study, the optimal operation modes of EM38 and EM31 were first determined to establish multiple linear regression models for estimating salinity from apparent soil electrical conductivity (ECa). Then, spatial trend and semi-variogram were analyzed, and spatial distribution of field salinity status was mapped. Their results suggested that ECa (EM38 and EM31) data correlated highly with salinity. Likewise, Ding and Yu [35] examined seasonal and spatial soil salinity changes using the electromagnetic induction method and field survey data in the Werigan-Kuqa River Delta Oasis, Xinjiang, China. Their results suggested that ECa obtained from EM38 was highly correlated with soil salinity.
Similarly, Hesham and Abdel-Rahman [36] used Geonics EM38 with field survey data to assess soil salinity in North Nile Delta, Egypt. In this study, statistical relations between ECa and ECe were first examined and then applied to predict and map salinity. They found that the actual and the predicted ECe levels were highly correlated. Besides, Cassel et al. [37] used a dual dipole EM38 device and the ground-truth EC survey data to quantify salinization's spatial and temporal progression and identify yield potential for tomatoes in California, USA. They found that the EM-predicted conductivity was consistent with the ground-truth EC survey data and increased with depth. Around 43% surveyed area was conducive to attaining 80% of the total yield potential of tomatoes.
With this research, we aimed to predict and map soil salinity severity in the subdistrict of the Non Thai district, Nakhon Ratchasima province, with a suitable interpolation method, which is a crucial method to estimate unknown data using known sample data [38]. Nielsen and Wendroth [39] categorized the interpolation method for digital soil mapping into four groups: non-geostatistical method, univariate geostatistical method, multivariate geostatistical method, and mixed method. In this study, the Inverse Distance At a national level, Khongnawang et al. [27] used an EM device, EM38, with field survey data and inversion software to map the distribution of clay and cation exchange capacity (CEC) in sandy and infertile soil in Ban Haet district, Khon Kaen, Thailand. In their study, they firstly examined a linear regression relationship between EM38 ECa with either topsoil (0-0.3 m), subsurface (0.3-0.6 m), or subsoil (0.6-0.9 m) clay and CEC. Then, the predicted equations were applied to map clay and CEC on three different levels.
Meanwhile, at an international level, Yao and Yang [34] used a mobile electromagnetic induction system, including Geonics EM38 and EM31, to perform an electromagnetic field survey and quantify the spatial pattern of soil salinity in Dongyihe Village of Kenli County, Shandong Province, China. In their study, the optimal operation modes of EM38 and EM31 were first determined to establish multiple linear regression models for estimating salinity from apparent soil electrical conductivity (ECa). Then, spatial trend and semi-variogram were analyzed, and spatial distribution of field salinity status was mapped. Their results suggested that ECa (EM38 and EM31) data correlated highly with salinity. Likewise, Ding and Yu [35] examined seasonal and spatial soil salinity changes using the electromagnetic induction method and field survey data in the Werigan-Kuqa River Delta Oasis, Xinjiang, China. Their results suggested that ECa obtained from EM38 was highly correlated with soil salinity.
Similarly, Hesham and Abdel-Rahman [36] used Geonics EM38 with field survey data to assess soil salinity in North Nile Delta, Egypt. In this study, statistical relations between ECa and ECe were first examined and then applied to predict and map salinity. They found that the actual and the predicted ECe levels were highly correlated. Besides, Cassel et al. [37] used a dual dipole EM38 device and the ground-truth EC survey data to quantify salinization's spatial and temporal progression and identify yield potential for tomatoes in California, USA. They found that the EM-predicted conductivity was consistent with the ground-truth EC survey data and increased with depth. Around 43% surveyed area was conducive to attaining 80% of the total yield potential of tomatoes.
With this research, we aimed to predict and map soil salinity severity in the subdistrict of the Non Thai district, Nakhon Ratchasima province, with a suitable interpolation method, which is a crucial method to estimate unknown data using known sample data [38]. Nielsen and Wendroth [39] categorized the interpolation method for digital soil mapping into four groups: non-geostatistical method, univariate geostatistical method, multivariate geostatistical method, and mixed method. In this study, the Inverse Distance Weighting (IDW), Ordinary Kriging (OK), Ordinary CoKriging (OCK), and Regression Kriging (RK) methods, which represent four categories, respectively, were selected to interpolate a continuous surface for soil salinity prediction.
Our specific research objectives were (1) to predict soil salinity using spatial interpolation methods and (2) to identify a suitable spatial interpolation method for soil salinity severity mapping. The proposed research framework can reduce the cost and time required for soil salinity surveying and mapping compared with the traditional methods of government agencies in Thailand. For instance, the LDD conducts a traditional survey by collecting soil profile information to map saline soil and its severity [14]. This traditional method requires a lot of time and resources. In addition, the distribution of soil salinity and its severity can be applied as primary information to solve the existing saline soil problem for agricultural uses.

Study Area
The study area is in the Non Thai district of Nakhon Ratchasima province, Thailand, between 15 • 4 18 N and 15 •  terpolate a continuous surface for soil salinity prediction.
Our specific research objectives were (1) to predict soil salinity using spatial interpolation methods and (2) to identify a suitable spatial interpolation method for soil salinity severity mapping. The proposed research framework can reduce the cost and time required for soil salinity surveying and mapping compared with the traditional methods of government agencies in Thailand. For instance, the LDD conducts a traditional survey by collecting soil profile information to map saline soil and its severity [14]. This traditional method requires a lot of time and resources. In addition, the distribution of soil salinity and its severity can be applied as primary information to solve the existing saline soil problem for agricultural uses.

Study Area
The study area is in the Non Thai district of Nakhon Ratchasima province, Thailand, between 15°4′18″ N and 15°18′34″ N, and 102°1′8″ E and 102°7′47″ E. This area is in the lower part of the Khorat Basin, supported by the rock salt and potash of the Mahasarakham formation. The Non Thai district covers ten sub-districts, namely, Non Thai, Dan Chak, Kampang, Samrong, Khang Phlu, Ban Wang, Banlang, Sai O, Thanon Pho, and Makha, with a total area of 547 sq. km ( Figure 2).
For over 30 years, the investors have operated salt mining and shrimp farming in more than 30 sites in the Non Thai, Non Sung, and Phra Thong Kham districts. Farmers have complained about some sites, but the investors did not comply with the Ministry of Industry's requirements. Furthermore, illegal salt mines have created many problems, mainly brine leakage into natural rivers and abandoned paddy fields.
Recently, the Nakhon Ratchasima Provincial Office [40] established a policy to solve the salinity problems occurring in the Non Thai district for arable land of farmers.

Materials and Methods
The research methodology entailed (1) apparent electrical conductivity measurement using an electromagnetic induction method, (2) in situ soil sample collection and electrical conductivity of the saturated soil paste extract measurement, (3) soil electrical conductivity estimation using linear regression analysis, (4) soil salinity prediction and accuracy assessment, (5) suitable interpolation method identification for soil salinity prediction, For over 30 years, the investors have operated salt mining and shrimp farming in more than 30 sites in the Non Thai, Non Sung, and Phra Thong Kham districts. Farmers have complained about some sites, but the investors did not comply with the Ministry of Industry's requirements. Furthermore, illegal salt mines have created many problems, mainly brine leakage into natural rivers and abandoned paddy fields.
Recently, the Nakhon Ratchasima Provincial Office [40] established a policy to solve the salinity problems occurring in the Non Thai district for arable land of farmers.

Materials and Methods
The research methodology entailed (1) apparent electrical conductivity measurement using an electromagnetic induction method, (2) in situ soil sample collection and electrical conductivity of the saturated soil paste extract measurement, (3) soil electrical conductiv-ity estimation using linear regression analysis, (4) soil salinity prediction and accuracy assessment, (5) suitable interpolation method identification for soil salinity prediction, and (6) soil salinity severity classification and overlay analysis with relevant data. See workflow in Figure 3. Details of each step are separately described in the following sections. Appl. Sci. 2022, 12, x FOR PEER REVIEW 5 of 35 and (6) soil salinity severity classification and overlay analysis with relevant data. See workflow in Figure 3. Details of each step are separately described in the following sections.

Apparent Electrical Conductivity Measurement Using an Electromagnetic Induction Method
Apparent electrical conductivity (ECa) data were measured using an EM38 sensor at systematic unaligned grid locations of 45 sites (Figure 4a). The number of sites was estimated as suggested by Kheoruenromne [41], who also recommends that the intensity of soil profile samples for a detailed reconnaissance soil survey at the scale of 1:40,000-1:100,000 should be one site per 2 sq. km because these scales and order of soil survey are fitted with the scale of DEM and geology data. The ECa was measured along three 100 m long transect lines at an interval of 50 m in nine locations at each site ( Figure 4b) during 19-21 August 2021. Thus, the total measured ECa for 45 sites is 405 (9 × 45) points. At each point, we used a dual dipole Geonics EM38 sensor, which provides simultaneous measurements with a horizontal measuring mode (HH) at a depth between 0 and 0.75 m and a vertical measuring mode (VV) at a depth between 0 and 1.5 m to measure ECa.

Apparent Electrical Conductivity Measurement Using an Electromagnetic Induction Method
Apparent electrical conductivity (ECa) data were measured using an EM38 sensor at systematic unaligned grid locations of 45 sites (Figure 4a). The number of sites was estimated as suggested by Kheoruenromne [41], who also recommends that the intensity of soil profile samples for a detailed reconnaissance soil survey at the scale of 1:40,000-1:100,000 should be one site per 2 sq. km because these scales and order of soil survey are fitted with the scale of DEM and geology data. The ECa was measured along three 100 m long transect lines at an interval of 50 m in nine locations at each site ( Figure 4b) during 19-21 August 2021. Thus, the total measured ECa for 45 sites is 405 (9 × 45) points. At each point, we used a dual dipole Geonics EM38 sensor, which provides simultaneous measurements with a horizontal measuring mode (HH) at a depth between 0 and 0.75 m and a vertical measuring mode (VV) at a depth between 0 and 1.5 m to measure ECa. Appl. Sci. 2022, 12,

In Situ Soil Sample Collection and Electrical Conductivity of the Saturated Soil Paste Extract Measurement
During the field survey, 30 samples from topsoil of about 0-25 cm depth and subsoil of about 25-50 cm depth at ten sites representing ten sub-districts were collected for soil electrical conductivity (ECe) extraction in the laboratory. In practice, soil samples at two levels at three locations of the middle transect line were collected using a handy auger. See Figure 4a,b. Later, fine soil samples were mixed with water at a ratio of 1:1 for measuring ECe with the electric conductivity meter based on the saturated soil paste extract method [42,43]. The measured ECe from the topsoil and subsoil layers were average for soil electrical conductivity estimation using linear regression analysis (LRA) as:

Soil Electrical Conductivity Estimation Using Linear Regression Analysis
Thirty average measured ECe data from the topsoil and subsoil samples and the observed ECa data using EM38 sensor at the exact locations (points) were first applied to identify a linear relationship using LRA. The performance of the LRA was examined to validate with three statistical values, correlation coefficient (R), coefficient of determination (R 2 ), and p-value, as suggested by [44][45][46] (see Table 2). Ramsey and Schafer [47] and Hurlbert and Lombardi [48] stated that the p-value is a continuous measure of the strength of evidence against the null hypothesis, with very small values indicating strong evidence of a difference between means (in a two-sample comparison), large values indicating little or no evidence of a difference.

In Situ Soil Sample Collection and Electrical Conductivity of the Saturated Soil Paste Extract Measurement
During the field survey, 30 samples from topsoil of about 0-25 cm depth and subsoil of about 25-50 cm depth at ten sites representing ten sub-districts were collected for soil electrical conductivity (ECe) extraction in the laboratory. In practice, soil samples at two levels at three locations of the middle transect line were collected using a handy auger. See Figure 4a,b. Later, fine soil samples were mixed with water at a ratio of 1:1 for measuring ECe with the electric conductivity meter based on the saturated soil paste extract method [42,43]. The measured ECe from the topsoil and subsoil layers were average for soil electrical conductivity estimation using linear regression analysis (LRA) as:

Soil Electrical Conductivity Estimation Using Linear Regression Analysis
Thirty average measured ECe data from the topsoil and subsoil samples and the observed ECa data using EM38 sensor at the exact locations (points) were first applied to identify a linear relationship using LRA. The performance of the LRA was examined to validate with three statistical values, correlation coefficient (R), coefficient of determination (R 2 ), and p-value, as suggested by [44][45][46] (see Table 2). Ramsey and Schafer [47] and Hurlbert and Lombardi [48] stated that the p-value is a continuous measure of the strength of evidence against the null hypothesis, with very small values indicating strong evidence of a difference between means (in a two-sample comparison), large values indicating little or no evidence of a difference.  [46] After that, the derived equations from HH and VV modes were separately used to estimate soil ECe data for 405 points. To avoid the negative values of ECe, the estimated ECe values of HH and VV modes were normalized using the following equation: Here, Corrected ECe is the ECe value of 405 points after correction, Estimated ECe is the estimated ECe value of 405 points from a simple linear equation, and Minimum value is the minimum value of 405 estimated ECe data. Later, the corrected 405 ECe data of HH and VV modes were randomly divided into the modeling dataset (70%), or 284 points, and the testing dataset (30%), or 121 points, for soil salinity prediction and accuracy assessment, respectively, in the next step.

Soil Salinity Prediction and Accuracy Assessment
The corrected ECe modeling dataset of HH and VV modes was used separately to interpolate a continuous surface for soil salinity prediction using IDW, OK, OCK, and RK methods.
Firstly, the IDW interpolation determines cell values using a linearly weighted combination of sample points. The weight is a function of inverse distance. The surface being interpolated should be that of a locationally dependent variable. This method assumes that the mapped variable decreases influence with distance from its sampled location [38]. The general equation of IDW by Shepard [49] is shown below.
Here,Ẑ(u 0 ) is the value being predicted for the target location u 0 ; N is the number of measured data points in the search window; w i represents the weights assigned to each measured point, and Z(u 1 ) is the observed value at location u i . u i = (x i ,y i ) [49]. Secondly, the OK method is a geostatistical interpolation method that weights the surrounding measured values to predict an unsampled location by incorporating the spatial dependence expressed by the semi-variogram in the estimation procedure [50]. The general equation of OK is displayed in the followings.
Here, Z(x 0 ) is the value of variable Z at the unsampled location x 0 ; λ i represents the weights associated with the data points, which consider the spatial relationship of the sampled points, and Z(x i ) is the observed value of Z at sampled locations x i [50]. Thirdly, the OCK method is a variation of the OK, which considers a secondary variable related to the primary variable that wants to estimate, in which additional observed variables (which are often correlated with each other and the variable of interest) are used to enhance the precision of the interpolation of the variable of interest at each location [51]. The general equation of OCK by McBratney and Webster [52] is presented in the following Equation (5).
Here, Z(x 0 ) is the position of the sample point, w 1i and w 2i are two regionalized variables, and z 1 x j and z 2 x j are weight coefficients [52]. In this study, the covariable factor of the OCK is soil moisture content applied for HH and VV modes. For soil moisture content measurement, some collected soil samples for ECe extraction in 30 locations were used to calculate soil moisture content based on a gravimetric basis. The soil moisture content (θ g ) is measured by weighting a soil sample (m wet ), drying the sample to remove the water, then weighting the dried soil (m dry ) [53] as: Figure 5 displays the spatial distribution of soil moisture content data, which was interpolated using IDW. The interpolated soil moisture content delivered a mean prediction error (ME) of 0.6403 percent by mass and root mean square error (RMSE) of 3.1008 percent by mass under the cross-validation of IDW. basis. The soil moisture content (θ ) is measured by weighting a soil sample (m ), drying the sample to remove the water, then weighting the dried soil (m ) [53] as: Figure 5 displays the spatial distribution of soil moisture content data, which was interpolated using IDW. The interpolated soil moisture content delivered a mean prediction error (ME) of 0.6403 percent by mass and root mean square error (RMSE) of 3.1008 percent by mass under the cross-validation of IDW. Fourthly, the RK method is a spatial interpolation technique that combines a regression of dependent variables on predictors with Kriging of the prediction residuals [54]. The general equation of RK by Odeh et al. [55] is shown in Equation (7).
Here, f(Q, x) is a function describing the structural component of S as a function of Q at x, and e (x) is the locally varying, spatially dependent residual from f (Q, x). In regression Kriging, the soil property S at an unvisited site is first predicted by f (Q, x), followed by Kriging of the model's residuals [55]. The RK method with a linear function and without a covariable factor was applied to predict ECe. In this study, IDW, OK, and OCK methods were conducted using ESRI ArcMap software, and the RK method was conducted with SAGA software.
After that, four standard statistical measurements, (1) mean prediction error (ME), (2) root mean square error (RMSE), (3) coefficient of determination (R 2 ), and (4) percent of bias (PBIAS), were used to assess the accuracy of each interpolation method based on the testing dataset, with 121 points, as suggested by Yao et al. [56] and Rafix et al. [57] Mathematically, the general form of ME, RMSE, R 2 , and PBIAS can be expressed in the following equations:  Fourthly, the RK method is a spatial interpolation technique that combines a regression of dependent variables on predictors with Kriging of the prediction residuals [54]. The general equation of RK by Odeh et al. [55] is shown in Equation (7).
Here, f(Q, x) is a function describing the structural component of S as a function of Q at x, and e (x) is the locally varying, spatially dependent residual from f (Q, x). In regression Kriging, the soil property S at an unvisited site is first predicted by f (Q, x), followed by Kriging of the model's residuals [55]. The RK method with a linear function and without a covariable factor was applied to predict ECe. In this study, IDW, OK, and OCK methods were conducted using ESRI ArcMap software, and the RK method was conducted with SAGA software.
After that, four standard statistical measurements, (1) mean prediction error (ME), (2) root mean square error (RMSE), (3) coefficient of determination (R 2 ), and (4) percent of bias (PBIAS), were used to assess the accuracy of each interpolation method based on the testing dataset, with 121 points, as suggested by Yao et al. [56] and Rafix et al. [57] Mathematically, the general form of ME, RMSE, R 2 , and PBIAS can be expressed in the following equations: Appl. Sci. 2022, 12, 10550 9 of 35 where n is the number of measured and predicted points.

Suitable Interpolation Method Identification for Soil Salinity Prediction
The derived four standard statistical measurements in the previous step, including ME, RMSE, R 2 , and PBIAS, were applied to compare and identify the suitable interpolation method for soil salinity prediction using the rank-sum technique of the Multi-Criteria Decision Analysis [58]. The spatial interpolation method that provides the highest total score was suitable for soil salinity prediction with estimated ECe data.

Soil Salinity Severity Classification and Overlay Analysis with Relevant Data
The interpolated soil salinity maps from suitable interpolation methods were further classified into five severity classes based on FAO standards, shown in Table 3. Moreover, overlay analysis was applied to describe the spatial relationship between soil salinity severity classification and relevant data, including geologic unit, soil texture, land use, and potential groundwater availability. Table 3. Classification of soil salinity severity based on FAO standards [59].

Soil Salinity Severity Class ECe (mS/m) Effect on Crop Plants
Non-saline 0-200 Salinity effects negligible Slightly saline 200-400 Yields of sensitive crops may be restricted Moderately saline 400-800 Yields of many crops are restricted Strongly saline 800-1600 Only tolerant crops yield satisfactorily Very strongly saline >1600 Only a few very tolerant crops yield satisfactorily Figure 6 shows an example of the measured ECa of the salt patch, including a field photograph and 2D and 3D maps from HH and VV modes by Surfer software, at Ban Don Nam Sai, Ban Wang sub-district.

Apparent Soil Electrical Conductivity Measurement
For HH mode, as representative of topsoil, the high ECa is found in the upper left corner in Figure 6b,c. Meanwhile, the high ECa is found in the western part in Figure 6d,e for VV mode, as representative of subsoil. This observation indicates the spatial variation of ECa's intensity in different soil depths.
Besides, Table 4 displays another example of the measured ECa with the essential statistical data of HH and VV modes in five land use types. Besides, Table 4 displays another example of the measured ECa with the essential statistical data of HH and VV modes in five land use types.    Table 4 shows the average measured ECa values of HH and VV modes in cassava and corn fields located in non-saline soil according to the FAO standard values. The average measured ECa values of bare land and paddy fields from the two modes were in slightly saline soil. On the contrary, average ECa values of the salt patch from the two modes were strongly saline.

In Situ Electrical Conductivity of the Saturated Soil Paste Extract Estimation
The laboratory results of the ECe from 30 plots are reported in Table 5. The average ECe values from two soil depth levels (0-25 and 25-50 cm) using Equation (1) vary from 15.84 mS/m at Sai O sub-district to 4359.5 mS/m at Kampang sub-district. The average ECe value of 30 plots is 1095.12 mS/m, and the standard deviation is 1291.09 mS/m. A large standard deviation value suggests that the measured values are scattered widely about the mean [60]. The average value of ECe by sub-district indicates strongly saline (800-1600 mS/m) and very strongly saline (>1600 mS/m) soil in six sub-districts, namely, Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang. This finding is reconfirmed by the spatial patterns of average ECe, which were interpolated based on average ECe data from 30 plots using the IDW method, as shown in Figure 7. Table 5. The saturated soil paste extraction (ECe) at the exact soil sample collection location.   Table 4 shows the average measured ECa values of HH and VV modes in cassava and corn fields located in non-saline soil according to the FAO standard values. The average measured ECa values of bare land and paddy fields from the two modes were in slightly saline soil. On the contrary, average ECa values of the salt patch from the two modes were strongly saline.

In Situ Electrical Conductivity of the Saturated Soil Paste Extract Estimation
The laboratory results of the ECe from 30 plots are reported in Table 5. The average ECe values from two soil depth levels (0-25 and 25-50 cm) using Equation (1) vary from 15.84 mS/m at Sai O sub-district to 4359.5 mS/m at Kampang sub-district. The average ECe value of 30 plots is 1095.12 mS/m, and the standard deviation is 1291.09 mS/m. A large standard deviation value suggests that the measured values are scattered widely about the mean [60]. The average value of ECe by sub-district indicates strongly saline (800-1600 mS/m) and very strongly saline (>1600 mS/m) soil in six sub-districts, namely, Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang. This finding is reconfirmed by the spatial patterns of average ECe, which were interpolated based on average ECe data from 30 plots using the IDW method, as shown in Figure 7.   Table 4 shows the average measured ECa values of HH and VV modes in cassava and corn fields located in non-saline soil according to the FAO standard values. The average measured ECa values of bare land and paddy fields from the two modes were in slightly saline soil. On the contrary, average ECa values of the salt patch from the two modes were strongly saline.

In Situ Electrical Conductivity of the Saturated Soil Paste Extract Estimation
The laboratory results of the ECe from 30 plots are reported in Table 5. The average ECe values from two soil depth levels (0-25 and 25-50 cm) using Equation (1) vary from 15.84 mS/m at Sai O sub-district to 4359.5 mS/m at Kampang sub-district. The average ECe value of 30 plots is 1095.12 mS/m, and the standard deviation is 1291.09 mS/m. A large standard deviation value suggests that the measured values are scattered widely about the mean [60]. The average value of ECe by sub-district indicates strongly saline (800-1600 mS/m) and very strongly saline (>1600 mS/m) soil in six sub-districts, namely, Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang. This finding is reconfirmed by the spatial patterns of average ECe, which were interpolated based on average ECe data from 30 plots using the IDW method, as shown in Figure 7.   Table 4 shows the average measured ECa values of HH and VV modes in cassava and corn fields located in non-saline soil according to the FAO standard values. The average measured ECa values of bare land and paddy fields from the two modes were in slightly saline soil. On the contrary, average ECa values of the salt patch from the two modes were strongly saline.

In Situ Electrical Conductivity of the Saturated Soil Paste Extract Estimation
The laboratory results of the ECe from 30 plots are reported in Table 5. The average ECe values from two soil depth levels (0-25 and 25-50 cm) using Equation (1) vary from 15.84 mS/m at Sai O sub-district to 4359.5 mS/m at Kampang sub-district. The average ECe value of 30 plots is 1095.12 mS/m, and the standard deviation is 1291.09 mS/m. A large standard deviation value suggests that the measured values are scattered widely about the mean [60]. The average value of ECe by sub-district indicates strongly saline (800-1600 mS/m) and very strongly saline (>1600 mS/m) soil in six sub-districts, namely, Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang. This finding is reconfirmed by the spatial patterns of average ECe, which were interpolated based on average ECe data from 30 plots using the IDW method, as shown in Figure 7.   Table 4 shows the average measured ECa values of HH and VV modes in cassava and corn fields located in non-saline soil according to the FAO standard values. The average measured ECa values of bare land and paddy fields from the two modes were in slightly saline soil. On the contrary, average ECa values of the salt patch from the two modes were strongly saline.

In Situ Electrical Conductivity of the Saturated Soil Paste Extract Estimation
The laboratory results of the ECe from 30 plots are reported in Table 5. The average ECe values from two soil depth levels (0-25 and 25-50 cm) using Equation (1) vary from 15.84 mS/m at Sai O sub-district to 4359.5 mS/m at Kampang sub-district. The average ECe value of 30 plots is 1095.12 mS/m, and the standard deviation is 1291.09 mS/m. A large standard deviation value suggests that the measured values are scattered widely about the mean [60]. The average value of ECe by sub-district indicates strongly saline (800-1600 mS/m) and very strongly saline (>1600 mS/m) soil in six sub-districts, namely, Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang. This finding is reconfirmed by the spatial patterns of average ECe, which were interpolated based on average ECe data from 30 plots using the IDW method, as shown in Figure 7.  Table 4 shows the average measured ECa values of HH and VV modes in cassava and corn fields located in non-saline soil according to the FAO standard values. The average measured ECa values of bare land and paddy fields from the two modes were in slightly saline soil. On the contrary, average ECa values of the salt patch from the two modes were strongly saline.

In Situ Electrical Conductivity of the Saturated Soil Paste Extract Estimation
The laboratory results of the ECe from 30 plots are reported in Table 5. The average ECe values from two soil depth levels (0-25 and 25-50 cm) using Equation (1) vary from 15.84 mS/m at Sai O sub-district to 4359.5 mS/m at Kampang sub-district. The average ECe value of 30 plots is 1095.12 mS/m, and the standard deviation is 1291.09 mS/m. A large standard deviation value suggests that the measured values are scattered widely about the mean [60]. The average value of ECe by sub-district indicates strongly saline (800-1600 mS/m) and very strongly saline (>1600 mS/m) soil in six sub-districts, namely, Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang. This finding is reconfirmed by the spatial patterns of average ECe, which were interpolated based on average ECe data from 30 plots using the IDW method, as shown in Figure 7.   Table 6 summarizes input data for the LRA to estimate the ECe of HH and VV modes. A scatter plot between the dependent variable (ECe average) and independent variables of two modes (ECa-HH and ECa-VV) for LRA is displayed in Figure 8. The results of LRA Figure 7. Spatial distribution of average ECe using the IDW method. Table 6 summarizes input data for the LRA to estimate the ECe of HH and VV modes. A scatter plot between the dependent variable (ECe average) and independent variables of two modes (ECa-HH and ECa-VV) for LRA is displayed in Figure 8. The results of LRA for soil ECe estimation of two modes are summarized in Table 7.    As a result, in Table 6 and Figure 8, it can be observed that ECa and ECe data can be divided into two groups. First, relatively large conductivities, where ECe is significantly greater than ECa, occur in the Non Thai, Dan Chak, Samrong, Kampang, Khang Phlu, and Ban Wang sub-districts. Second, relatively small conductivities, where ECa is significantly greater than ECe, occur in the Sai O, Makha, Thanon Pho, and Banlang sub-districts. The main reason to support these phenomena is the variation of soil conductivities due to the dynamic nature of the soluble salts in the study area.

Linear Regression Analysis for Soil Electrical Conductivity Estimation
Moreover, it can be observed that 13 samples with ECe values less than 100 show a moderate correlation with ECa values. These samples should be considered separately. However, when we separate data into two groups, the number of samples is too small. Jenkins and Quintana-Ascencio [61] suggested that the minimum sample size for linear regression should be equal to or more than 25. Additionally, this study aims to examine the distribution of soil salinity and its severity in the whole study area.
As resulted in Table 7, the R value between ECa and ECe of the two modes are 0.873 and 0.861, respectively. ECa has a very strong positive relationship with ECe [44,62]. The R 2 values of HH and VV modes are 0.762 and 0.742, respectively. Both R 2 values provide a very strong goodness-of-fit measure for the linear regression model [45,63,64]. Besides, the p-values of regression of HH and VV are 0.000 and 0.000, respectively. The derived p-values of the two modes, with a value of less than 0.01, show very strong evidence of a difference between means in a two-sample comparison [46,65,66]. So, the result of LRA can be validated and accepted for estimating ECe in the remaining ECa locations.

Soil Salinity Prediction and Accuracy Assessment
The results of soil salinity prediction and accuracy assessment using four spatial interpolation methods (IDW, OK, OCK, and RK) of two measuring modes are separately summarized in the following sections.

Soil Salinity Prediction and Accuracy Assessment Using the IDW Method
The spatial distribution maps of soil salinity prediction of HH and VV modes using IDW with a power parameter of 2 are displayed in Figure 9. Essential statistical data of soil salinity prediction with ECe data using IDW of two modes by sub-district are reported in Table 8. soil salinity prediction with ECe data using IDW of two modes by sub-district are reported in Table 8.   In Figure 9, the patterns of soil salinity distribution from two modes at the sub-district levels are slightly different. Under the HH mode, at a depth of 0-0. 75  In addition, the result of the accuracy assessment of the IDW method with the testing dataset was reported in Table 9. As a result, the ME and RMSE values of the HH mode are about 663.18 and 1308.51 mS/m, and those of the VV mode are about −80.63 and 803.92 mS/m. The R 2 of the HH mode is 0.8422, and that of the VV mode is 0.8354. The PBIAS of the HH mode is about −43.67%, and that of the VV mode is about 3.23%. Table 9. Statistical data for accuracy assessment of the IDW method. When ME values are compared between two modes, the HH mode delivers a higher error in the mean. Likewise, the RMSE value of the HH mode is higher than the VV mode. Moreover, the PBIAS value of HH mode is higher than VV mode. In the meantime, the R 2 values of HH and VV modes are indifferent, providing a very strong correlation between measured and predicted values.

Soil Salinity Prediction and Accuracy Assessment Using the OK Method
The spatial distribution maps of soil salinity prediction of two modes using the OK method with defined parameters on semi-variogram modeling and standard searching neighborhood are displayed in Figure 10. Essential statistical data of soil salinity prediction using OK of two modes by sub-district are reported in Table 10. The spatial distribution maps of soil salinity prediction of two modes using the OK method with defined parameters on semi-variogram modeling and standard searching neighborhood are displayed in Figure 10. Essential statistical data of soil salinity prediction using OK of two modes by sub-district are reported in Table 10.   As shown in Figure 10, the variations in soil salinity distribution from two measuring levels of two modes at the sub-district level are slightly different. In Table 10 In addition, the result of the accuracy assessment of the OK method with the testing dataset was reported in Table 11. As a result, the ME and RMSE values of the HH mode are about 632.85 and 1247.72 mS/m, and those of the VV mode are about −145.42 and 807.13 mS/m. The R 2 value of the HH mode is 0.8526, and that of the VV mode is 0.8344. The PBIAS of the HH mode is about −41.68%, and that of the VV mode is about 5.83%. When comparing ME values between two measuring modes, the VV has a lower mean error. Likewise, the RMSE and PBIAS values of VV are lower than HH. In the meantime, the R 2 values of HH and VV modes are slightly different, providing a very strong correlation between measured and predicted values.

Soil Salinity Prediction and Accuracy Assessment Using OCK Method
Spatial distribution maps of soil salinity prediction of two measuring modes using OCK with covariable factor, soil moisture content, and the defined parameters on semivariogram modeling and standard searching neighborhood are displayed in Figure 11.
Primary statistical data of soil salinity prediction using OCK of two measuring modes by sub-district are reported in Table 12.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 19 of 35 Primary statistical data of soil salinity prediction using OCK of two measuring modes by sub-district are reported in Table 12.
(a) (b) Figure 11. Spatial distribution of soil salinity prediction using OCK of two measuring modes: (a) HH and (b) VV.  As shown in Figure 11, patterns of soil salinity distribution from the two modes at the sub-district level are slightly different. At the same time, in Table 12  When comparing ME values between two measuring modes, VV has a lower mean error. Similarly, the RMSE value of VV is lower than HH. Likewise, the PBIAS value of VV has a lower magnitude. In the meantime, the R 2 values of HH and VV modes are slightly different, providing a very strong correlation between measured and predicted values.

Soil Salinity Prediction and Accuracy Assessment Using the RK Method
The spatial distribution maps of soil salinity prediction of two measuring modes using RK with stepwise regression method and ordinary Kriging type are displayed in Figure 12. Essential statistical data of soil salinity prediction using RK of two measuring modes by sub-district are reported in Table 14.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 21 of 35 Figure 12. Essential statistical data of soil salinity prediction using RK of two measuring modes by sub-district are reported in Table 14.
(a) (b)    In Figure 12, the patterns of soil salinity distribution from two measuring modes at the sub-district level are slightly different. In Table 14 Moreover, the results of the accuracy assessment of the RK method with the testing dataset are reported in Table 15. The ME and RMSE values of the HH mode are about 570.63 and 1318.15 mS/m, and those of the VV mode are about −420.68 and 1218.47 mS/m. The R 2 values of the HH and VV modes are 0.8280 and 0.6665, respectively. The PBIAS of the HH mode is about −37.58%, and that of the VV mode is 16.88%. When comparing ME values of two modes, VV has the lower error in the mean. Similarly, the RMSE value of VV is lower than that of HH. Likewise, the PBIAS value of VV is lower than that of HH. However, the R 2 value of VV is lower than that of HH and provides a moderate correlation between measured and predicted values.

Suitable Interpolation Method for Soil Salinity Prediction and Mapping
The comparison of the accuracy of four interpolated methods in terms of model performance for soil salinity prediction in the two measuring modes, including ME, RMSE, R 2 , and PBIAS, is reported in Table 16. As a result of ranking, the higher efficiency value yields a higher score for total score calculation using the rank-sum method, as reported in Table 17.  Table 16 shows the lowest ME value of the four methods for HH mode is OK and that for VV mode is IDW. The lowest RMSE value of the four methods for HH mode is the OCK, and that for the VV mode is IDW. The highest R 2 value of the four methods for HH mode is the OCK, and that for the VV mode is IDW. Besides, the lowest PBIAS value of the four methods for HH mode is OK, and that for the VV mode is IDW.
After applying the rank-sum method as a summary in Table 17, it can be seen that the OK and OCK methods deliver the highest total score, with 12. However, OCK provides the RMSE and R 2 values better than OK, while there is no difference in the ME and PBIAS values between OK and OCK. So, OCK is chosen as a suitable spatial interpolation method for soil salinity prediction of HH measuring mode. Meanwhile, the IDW is a suitable spatial interpolation method for soil salinity prediction of the VV measuring mode because it delivers the highest total score, with a value of 16. The derived maps of suitable spatial interpolation methods of the two measuring modes are further used to classify soil salinity severity classes.

Soil Salinity Severity Classification
The derived soil salinity prediction of the HH and VV modes using OCK and IDW in Figures 11a and 9b, respectively, were classified as soil salinity severity according to FAO standards (see Table 2), as shown in Figure 13. The area of soil salinity severity classification of the two measuring modes by sub-district and district is summarized in Table 18.    For the HH mode using OCK, the area of the non-saline class is only found in two sub-districts, Dan Chak and Ban Wang, and it covers an area of about 11.25 sq. km. The top three sub-districts with slightly to moderately saline classes are Non Thai, Banlang, and Dan Chak, with areas varying from 12.05 to 24.21 sq. km. In contrast, three sub-districts, Khang Phlu, Kampang, and Samrong, display the minor area of slightly to moderately saline classes varying from 1.59 to 8.25 sq. km. Meanwhile, the top three sub-districts with strongly to very strongly saline classes are Non Thai, Banlang, and Khang Phlu, with areas varying from 56.50 to 82.97 sq. km. On the contrary, three sub-districts, namely Dan Chak, Sai O, and Thanon Pho, display the minor area of strongly to very strongly saline classes, varying from 12.78 to 33.16 sq. km.
For the VV mode using OK, the area of the non-saline class is only found in one subdistrict, Dan Chak, and it covers an area of about 0.54 sq. km. The top three sub-districts with slightly to moderately saline classes are Makha, Ban Lang, and Dan Chak, with areas varying from 9.25 to 21.55 sq. km. In contrast, three sub-districts, namely, Samrong, Khang Phlu, and Sai O, display the minor area of slightly to moderately saline classes varying from 0.17 to 2.53 sq. km. Meanwhile, the top three sub-districts with strongly to very strongly saline classes are Non Thai, Banlang, and Khang Phlu, with areas varying from 62.36 to 83.80 sq. km. On the contrary, three sub-districts, Dan Chak, Sai O, and Thanon Pho, display the minor area of strongly to very strongly saline classes, varying from 26.14 to 33.16 sq. km.

Soil Salinity Severity Classification and Geologic Unit
The overlay analysis between soil salinity severity classification of two measuring modes and the geologic unit map of DMR in 2007 [67] is reported in Table 19. In Table 19, it can be seen that the most dominant soil salinity severity classification of two measuring modes over two main geologic units, the Maha Sarakham formation and Quaternary sediments, are strongly and very strongly saline classes. They cover areas in the HH and VV modes, about 81.48% and 89.62%, respectively.

Soil Salinity Severity Classification and Soil Texture
Based on soil series data of LDD in 2011 [68], soil textures and other areas in the study area are composed of loamy sand, sandy loam, silt loam, silty clay, and clay with a salt farm, urban areas, and water bodies. The result of overlay analysis between soil salinity severity classification of two measuring modes and soil texture is reported in Table 20. In Table 20, for the result of the HH measuring mode, the top three dominant soil textures in the slightly saline and moderately saline classes are sandy loam, silty clay, and clay; they cover about 29.27, 25.14, and 11.81 sq. km, respectively. Meanwhile, the top three dominant soil textures in the strongly saline and very strongly saline classes are sandy loam, silty clay, and loamy sand, which cover about 254.79, 66.04, and 41.77 sq. km, respectively.
For the result of the VV measuring mode, the top three dominant soil textures in the slightly saline and moderately saline classes are sandy loam, silty clay, and silt loam. They cover about 20.94, 16.95, and 7.21 sq. km, respectively. In contrast, the top three dominant soil textures in the strongly saline and very strongly saline classes are sandy loam, silty clay, and loamy sand. They cover 271.03, 74.23, and 44.89 sq. km, respectively.

Soil Salinity Severity Classification and Land Use
The overlay analysis between soil salinity severity classification of two measuring modes and land use data of LDD in 2015 [69] is reported in Table 21. From the HH measuring mode data in Table 21, the top three dominant land use types in the slightly saline and moderately saline classes are paddy fields, field crops, and urban and built-up areas. They cover about 120.97, 46.48, and 12.90 sq. km, respectively. Likewise, the top three dominant land use types in the strongly saline and very strongly saline classes are paddy fields, field crops, and urban and built-up areas. They cover about 226.57, 58.70, and 17.63 sq. km, respectively. Additionally, soil salinity severity classification of the HH measuring mode in salt flats consists of moderately saline, strongly saline, and very strongly saline.
Based on the VV measuring mode, the top three dominant land use types in the slightly saline and moderately saline classes are paddy fields, field crops, and miscellaneous land. They cover about 33.77, 11.30, and 4.02 sq. km, respectively. Similarly, the top three dominant land use types in strongly saline and very strongly saline classes are paddy fields, field crops, and built-up areas. They cover about 326.17, 94.23, and 28.14 sq. km, respectively. Soil salinity severity classifications of the VV measuring mode in salt flats are mainly moderately saline, strongly saline, and very strongly saline. Table 22 contains the overlay analysis between soil salinity severity classification of two measuring modes and potential groundwater availability data of the Department of Groundwater Resources (DGR) in 2015 [70], consisting of the total dissolved solids (TDS) > 1500 mg/L with the expected well yield of <2 cu. m/h and TDS > 1500 mg/L with 2-10 cu. m/h.  Table 22 shows that the most dominant soil salinity severity classifications in the two measuring modes are strongly and very strongly saline classes. The areas of both classes in HH and VV modes are 81.76% and 89.70% of the study area. Meanwhile, slightly saline and moderately saline classes in the two modes are 16.50% and 10.30% of the study area, respectively.

Estimation of ECe of HH Mode and VV Mode Using Linear Regression Analysis
The ECe data of HH and VV modes, as the primary input for soil salinity prediction, were estimated based on the derived equation from linear regression analysis. The derived equations of HH and VV modes provided the R values, with a value of 0.873 and 0.861, the R 2 , 0.762 and 0.742, and the p-value, with a value of 0.000 and 0.000, respectively. As a result, ECa has a very strong positive relationship with ECe, and both R 2 can explain the variance in observations, ECe, about 76% and 74% [62]. The derived p-values of the two modes, with a value of less than 0.01, show very strong evidence of a difference between means in ECa and ECe sample comparison [47,48].
This study agrees with many similar studies [35,36,71]. Hesham and Abdel-Rahman [36] examined the statistical relationship between ECa and ECe for salinity predicting and mapping. They found that ECa values were highly correlated with ECe values, with R 2 > 0.70. Likewise, Jiang et al. [71] used an EM38 sensor to measure ECa in four land-use types: natural desert, natural vegetation, apple orchard, and winter wheat farmland, and to identify a linear relationship with 60 soil ECe samples. They found a strong linear relationship between observed ECa and measured ECe, with an R 2 value of 0.87. They used the derived equation to convert all observed ECa data into ECe data for spatial interpolation with the OK method to classify soil salinity severity.

Soil Salinity Prediction Using Interpolation Method
Salts are highly soluble in surface and groundwater and can be transported with water movement. There are two types of salinity: primary and secondary. Primary salinity occurs naturally in soils and waters, such as salt lakes, pans, marshes, and salt flats. Meanwhile, secondary salinity is salting from human activities, usually land development and agriculture [72] (https://omexcanada.com, accessed on 10 August 2022). Generally, soil salinity can vary temporally and spatially due to the dynamic nature of the soluble salts [73][74][75].
The primary and secondary types of salinity are found in the Non Thai district since it is situated over deep-seated rock salt [6,8] and has human activities. Based on the modeling dataset, the soil salinity prediction of two measuring modes was successfully interpolated using the IDW, OK, OCK, and RK methods. The spatial pattern of soil salinity distribution from two measuring modes with four different methods at the sub-district level was slightly different, see  The mean values of ECe using the four methods of the two measuring modes at the district level varied from about 2156 to 2293 mS/m for HH and from about 2377 to 2401 mS/m for VV. This result indicates a significant difference among soil spatial salinity variances from two measured soil depths (topsoil and subsoil). When comparing CV values between topsoil (HH node) and subsoil (VV mode) of four interpolation methods, the topsoil delivered a higher CV than that of the subsoil. See Tables 8, 10, 12 and 14. This result indicates a variation of soil salinity at the topsoil level that is higher than that at the subsoil level since most topsoil is utilized for agricultural uses, i.e., paddy fields and field crops. According to the 2015 Land Development Department report, paddy fields and field crops cover about 465.80 sq. km or about 88.22% of the study area [68]. The comprehensive effects of human activities are critical factors for the salt distribution process [76].

Accuracy Assessment of Interpolation Method
When comparing ME values of two modes from four interpolation methods, VV, representative of subsoil level, has a lower error in the mean than that HH, representative of topsoil. Similarly, the RMSE value of VV is lower than that of HH. Likewise, the PBIAS value of VV is lower than that of HH. However, there is no significant difference between the R 2 values of HH and VV. These results indicate that the predicted ECe data at the subsoil layer from the four methods deliver more accurate results than at the topsoil layer. Interestingly, all the predicted ECe data of the topsoil level of the four methods are overestimated. Still, all the predicted ECe data of the subsoil level of the four methods are underestimated.

Suitable Interpolation Method for Soil Salinity Prediction
The standard measurement of accuracy assessment (ME, RMSE, R 2 , and PBIAS), as suggested by [56,57] with the rank-sum technique [58], was applied to identify a suitable interpolation method for soil salinity prediction. As a result, OCK, as a multivariate geostatistical method [39], was suitable for soil salinity prediction at the topsoil level (HH mode). In a similar study, Abdennour et al. [28] conducted predictive mapping of soil electrical conductivity as a proxy of soil salinity in southeast Algeria by measuring 42 ECe samples collected from topsoil (0-15 cm). Three models were compared: OK, OCK, and IK (Indicator Kriging). They found that the OCK was the best, with the lowest ME and RMSE values, where the ME was 9 mS/m, and the RMSE was 153 mS/m. Meanwhile, a suitable method for soil salinity prediction at subsoil level (VV mode) was IDW as a non-geostatistical method [39]. This finding agrees with Hammam and Mohamed [77], who mapped soil salinity in the Eastern Nile Delta, Egypt, by analyzing estimated 92-point ECe values using cross-validation of IDW with power 3. The ME value was 12 mS/m, and the RMSE value was 39 mS/m. Nevertheless, the derived ME and RMSE values of OCK and IDW in this study were relatively high because the ECe data for soil salinity interpolation was not directly measured in the field. Still, they were estimated based on the derived equation of each measuring mode using LRA.

Soil Salinity Severity Classification
The spatial distribution of soil salinity severity classification using the suitable method in the two measuring modes differs depending on the estimated ECe data from the modeling dataset and the applied method. At the district level, the most dominant soil salinity severity classifications of the HH and VV modes are strongly and very strongly saline classes, with an area of about 81.50% and 89.62% of the whole area, respectively. This finding is expected since researchers and government agencies have recognized the saline soil problem in the Non Thai district [14,78].

Soil Salinity and Geologic Unit
The strongly and very strongly saline soil classes of HH and VV modes covered areas of two geologic units, the Maha Sarakham formation and Quaternary sediments, about 81% and 90%, respectively. This finding is consistent with the previous report about soil salinity in the Northeast region of Thailand [6][7][8][9][10][11]. Based on annual surveys from 1973 to 1982 by the Economic Geology Division, Department of Mineral Resources, the Maha Sarakham formation comprises evaporites and redbeds, attaining a probable original thickness of about 300 m. The rock units (member status) begin with the Basal Anhydrite, Lower Salt, Lower Clastic, Middle Salt, Middle Clastic, Upper Salt, and Upper Clastic at the top. Sedimentary features preserved in cores also point toward a nonmarine origin. This interpretation leads to inferences about the formation of saline soils in northeast Thailand in many cases, such as the continuity of the potash horizons and the occurrence and migration of subsurface brine [79].
Likewise, Wongsomsak [78] studied salinization in northeast Thailand using various exposures and analytical data from the Maha Sarakham formation's rock salt member at three points in the Non Thai district: (1) dug pond at Ban Khok Krasang, Non Thai district, in which were found red siltstone at a depth of 2.8 m with pH 9.5, ECe 24,000 mS/m, and Na total base of 6.5%; (2) slightly undulating rice field with termite mounds, 1 km west of Ban Khang Phlu, on the Khok Kruat-Non Thai road, in which were found white sandy clay, many calcium carbonate nodules at a depth of 1.2 m with pH 8.8, ECe 128,000 mS/m, and Na total base of 28.64%; (3) very slightly undulating rice field, north of Khok Sung village along with Cho Ho-Non Thai road, in which were found light olive gray, loamy clay at a depth of 2.2 m with pH 9.2, ECe 62,000 mS/m, and Na total base of 18.97%. All three had ECe values greater than 1600 mS/m and soil salinity severity classes of very strongly saline.

Soil Salinity and Soil Texture
The spatial relationship between soil salinity and texture is consistent with a report by Boonsompopphan et al. [80]. They studied how to identify soil series in the field, enabling Thai farmers to share their knowledge and expertise. Most of Thailand's inland salt-affected soils have a leached surface horizon of sandy loam or loamy sand overlying a very hard and impermeable Bt horizon (natric horizon). The subsoils are sandy clay loam or clay loam and are generally characterized by columnar or prismatic structure. Due to flat topography coupled with impervious layers, the soils show dominant signs of wetness. A gray color matrix with brownish or yellowish mottles is present throughout the soil profile. Salt and lime concretions, if present, are found in the subsoil. However, in the dry season, salt crusts are observable on the soil surface. In the profile, the dominant salt contains sodium, and its amount commonly increases with increasing depth.

Soil Salinity and Land Use
Agricultural practices play a critical role in salt salinity distribution [76]. The spatial relationship between soil salinity and land use types in this study indicates that paddy fields and field crops dominate the study area. This finding is similar to the previous studies. Madyaka [81] studied the spatial modeling and prediction of soil salinization using SaltMod in a GIS environment in the Nong Suang sub-district, Kham Thale So district, Nakhon Ratchasima, Thailand, by measuring ECe and pH in 51 study areas at three levels: 0-30, 30-60, and 60-100 cm. The SaltMod was used to model temporal changes in salinization over two decadal (20 years) periods. The pattern indicates the effect of physiographic conditions on salinity, as a major part of the southwestern side is dominated by ridges and the northeastern side by flood plains, lateral vales, and terraces. This pattern can also be associated with land use types, as the latter side is dominated by paddy fields, while cassava and maize are cultivated southwestern.
Hall et al. [82] studied land use and hydrological management with the Isaan Catchment Hydrogeological and Agricultural Model (ICHAM), an integrated model at a regional scale in Northeast Thailand. The ICHAM consists of a series of interconnected worksheets, which cover different aspects of the salinity management system. It brings knowledge of hydrogeology, agronomy, and farm economics, but further ground studies are needed before making confident recommendations about changing land use in particular areas. This finding shows how salinization can be managed, the general direction to take, and the approximate magnitude of land use change needed. The simulation results show that salinity will almost double in the next 30 years if current land use is maintained for a catchment in northeast Thailand.

Soil Salinity and Potential Groundwater Availability
The strongly and very strongly saline soil classes of two measuring modes covered areas of two potential groundwater availability classes, TDS) > 1500 mg/L with the expected well yield of <2 cu. m/h and TDS > 1500 mg/L with 2-10 cu. m/h, about 82% and 90%, respectively. This finding is consistent with a previous report about soil salinity in the northeast region of Thailand. According to the Bureau of Groundwater Conservation and Restoration, the groundwater supply region in the Non Thai district ranges from brackish to salt water. It has a high chloride (Cl) chemical composition and high total soluble TDS [70].
Likewise, Arunin [83] studied the reforestation of landscape salinity project conducted in the Nakhon Ratchasima province's potential salt source recharge area using an EM survey to identify and map the recharge and discharge areas and migration of subsurface brine. It was found that groundwater levels were low in the recharge area, ranging from 4.81 to 5.87 m below the soil surface, while in the discharge area, the groundwater levels were high, ranging from 0.48 to 1.04 m. The lower groundwater levels are attributed to the high consumptive water use of trees, which was supported by another finding using sap flux density and measurements of heat pulse velocity, finding that the annual water use of eucalyptus plantations amounted to 1230 and 270 mm/year in moderately and strongly saline soils, respectively.

Conclusions
Handheld electromagnetic sensors, EM38, were used to measure ECa at 405 locations with two standard measuring modes: HH (0-0.75 m depth) and VV (0-1.5 m depth) over 45 sites in the study area. Thirty soil samples from the topsoil and subsoil layers were collected for the electrical conductivity of the saturated soil paste extract measurement in the laboratory. Then, a linear relationship between the measured ECa and ECe data at the exact location was identified using linear regression analysis for estimating the remaining ECa data. The results showed a strong positive relationship between ECa and ECe. R values of HH and VV modes were 0.873 and 0.861, and R 2 were 0.762 and 0.742, respectively. The model has a p-value of less than 0.001. The total calibrated soil conductivity data were further randomly categorized into two datasets: the modeling dataset (70%), with 290 points, and the testing dataset (30%), with 123 points, for soil salinity prediction and validation, respectively.
For soil salinity prediction, the modeling dataset with 284 points was successfully applied to predict soil salinity of two measuring modes using IDW, OK, OCK, and RK methods. The patterns of soil salinity distribution from two measuring modes with four different methods at the sub-district level were slightly different. After that, the predicted results of the four methods were assessed for accuracy based on the testing dataset with 121 points using ME, RMSE, R 2 and PBIAS, and they were applied to identify a suitable interpolation method for soil salinity prediction. It was found that OCK is a suitable spatial interpolation method for soil salinity prediction in HH mode. At the same time, IDW is a suitable spatial interpolation method for soil salinity prediction for the VV mode.
Finally, as a suitable spatial interpolation method, the derived soil salinity prediction data of HH and VV modes using OCK and OK were further applied to classify soil salinity severity according to FAO standards. As expected, the most dominant soil salinity severity classifications under the two measuring modes were strongly and very strongly saline. The spatial distribution of soil salinity severity would be beneficial in supporting site-specific strategies for crop management. Consequently, this research methodology can be applied to conduct soil salinity surveys in other areas to reduce costs and time for government agencies in Thailand.
Nevertheless, to apply the EMI method for soil salinity surveys, the users should first understand the principle of EMI and how to calibrate and operate the EM device properly for accurate ECa measurement. Additionally, sophisticated regression analysis, such as non-linear regression, should be examined to identify a relationship between ECa and ECe for estimating the ECe of the remaining sample points.