Intercropping Simulation Using the SWAP Model: Development of a 2×1D Algorithm

Abstract

Intercropping is a common cultivation system in sustainable agriculture, allowing crop diversity and better soil surface exploitation. Simulation of intercropped plants with integrated soil–plant–atmosphere models is a challenging procedure due to the requirement of a second spatial dimension for calculating the soil water lateral flux. Evaluations of more straightforward approaches for intercrop modeling are, therefore, mandatory. An adaptation of the 1D model Soil, Water, Atmosphere and Plant coupled to the World Food Studies (SWAP/WOFOST) to simulate intercropping (SWAP 2×1D) based on solar radiation and water partitioning between plant strips was developed and the outcomes are presented. An application of SWAP 2×1D to maize–soybean (MS) strip intercropping was evaluated against the monocropping maize (M) and soybean (S) simulated with the 1D model SWAP/WOFOST, and a sensitivity analysis of SWAP 2×1D was carried out for the intercropping MS. SWAP 2×1D was able to simulate the radiation interception by both crops in the intercropping MS and also to determine the effect of the radiation attenuation by maize on soybean plants. Intercropped plants presented higher transpiration and resulted in lower soil evaporation when compared to their equivalent monocropping cultivation. A numerical issue involving model instability caused by the simulated lateral water flux in the soil from one strip to the other was solved. The most sensitive plant parameters were those related to the taller plant strips in the intercropping, and soil retention curve parameters were overall all significantly sensitive for the water balance simulation. This implementation of the SWAP model presents an opportunity to simulate strip intercropping with a limited number of parameters, including the partitioning of radiation by a well-validated radiation sharing model and of soil water by simulating the lateral soil water fluxes between strips in the 2×1D environment.

1. Introduction

Intercropping is a method employed for agricultural sustainability in which two or more plant species are simultaneously cultivated together on the same land [1,2,3] Intercropping systems are favored by a better land occupation, a decrease in the use of agrochemicals, better soil conservation, and yield diversification, among other benefits [4,5,6,7]. Management challenges for intercropping include the use of an optimum plant density, cultivation design or plant strip widths, and the best times of seeding for each of the crops to maximize crop productivities [8]. A better strategy of intercropping cultivation can be facilitated and accomplished with multi-species cropping system modeling [9].

Models for intercropping systems have been developed using several approaches [10] simulating the competition of intercrop plants mainly for water and radiation. The competition for radiation by intercropped plants is, however, the main process taken into account in modeling studies [11,12,13,14,15,16,17]. Chimonyo [18] pointed out the difficulties of evaluating the root systems and water extraction dynamics of multicropping systems are the main reason for the modeling of the processes belowground to be put aside. Another reason for avoiding the belowground modeling in intercropping is the difficulty of integrating two spatial dimensions for the simulation of the vertical and horizontal soil water fluxes [19].

Approaches in which the belowground resource sharing can be simulated simply and feasibly are highly desirable [20,21]. Chen [22] evaluated the concept that the lateral movement of water within strips of an intercropping occurs predominantly when a gradient of soil water content exists between strips. This procedure, physically based on Darcy’s law, was used also in Karray [23] for modeling the water exchange between plots of olive-trees and annual crops. Such a simple concept can facilitate the modeling of the lateral soil water flow in the strip intercropping with 1D-models in which only the vertical coordinate is simulated.

The hydrological model Soil, Water, Atmosphere and Plant (SWAP) [24] coupled to the crop growth model World Food Studies (WOFOST) [25] simulates the hydrological processes and plant growth in cropping systems and provides detailed water balances and crop productivity predictions. SWAP/WOFOST is a 1D model, simulating flow processes in the vertical direction. Results may be considered representative of a plot or a field-scale region if soil, cropping system, and meteorological data are considered spatially uniform. Implicitly, SWAP/WOFOST is unable to simulate two crops growing together and partitioning water and radiation, like in intercropping or agroforestry. Intercropping simulation with SWAP could be accomplished, however, by simulating two soil-plant modules simultaneously, having a third routine controlling the interchange of water and radiation between plants.

In this context, we aimed to develop a modified version of SWAP/WOFOST that allows simulating upper and belowground resources partitioning between intercropped plants. In the present evaluation, we considered that in each intercropping strip the plant and root system are limited to the respective soil profile as well as to the aerial portion of the strip. Hence, plants can be simulated separately in modules, sharing radiation in accordance to plant strip geometry, and interchanging water between modules by a lateral flow dependent on the soil hydraulic status of the soil layers in both strips. The new version of SWAP/WOFOST combines the robustness of the existing 1D model to the possibility of simulating the water balance components and crop growth in an intercropping scenario. Model theory and assumptions are described and tested for a hypothetical case of relay-strip intercropping.

2. Material and Methods

2.1. Model Description

2.1.1. SWAP 2×1D

The intercropping simulation model SWAP 2×1D consists of two joined conventional SWAP/WOFOST (1D) (version 4.0.4) modules, capable of interchanging data of plant and soil-water fluxes after every 24 h simulation period. The daily loops of simulation (the LOOP process in Figure 1) for the desired period of intercropping cultivation with SWAP 2×1D are executed by an external routine named Share in combination with a batch file. Reading and storage of data (INPUT process) takes place once at the onset of a simulation. Output data are written to default SWAP document files every loop until the end of the simulation.

2.1.2. Radiation Interception

A radiation interception model for relay-strip intercropping based on a canopy formed by block-structures was developed by Pronk [15], and evaluated by Gou [11,12]. The radiation interception by strip intercropping plants as described by Pronk [15] and Gou [11,12] remains in between the radiation intercepted by a homogeneous and a compressed canopy. In a homogeneous canopy (Figure 2), leaf areas are considered homogeneously distributed horizontally, like in mixed pastures. The radiation interception of a homogeneous canopy is higher than the radiation interception by a strip intercropping canopy for two main reasons: (1) A homogenous canopy model does not take into account the soil paths present in the relay-strip intercropping when only one crop species is in the field; and (2) the border effects of the strips are not taken into account in the homogenous canopy model. In a compressed canopy scenario, the plant rows and the soil paths can be thought of as two separate units: (1) plant rows without paths and (2) uncovered soil paths, both occupying a complementary fraction of the total soil area (Figure 2). In the compressed type canopy model, the radiation is partly intercepted by the canopy with plants very close to each other, and partly by the uncovered soil. The strip-canopy radiation interception is larger than in a compressed canopy because the strip-border plants are more exposed to the radiation and are more numerous than in the compressed canopy [11,12].

In a relay-strip intercropping, one of the canopies is taller than the other and referring to radiation interception the taller one is subdivided in an upper part (above the shorter canopy) and a lower part, according to the computation algorithm described in the following.

The fraction of the intercepted radiation by a strip-canopy (f_t) is in between a maximum value, the radiation interception fraction by a homogeneous canopy, and a minimum value equal to the addition interception fraction by a compressed canopy. f_t is calculated as the weighted average of the radiation interception fraction by a homogeneous and a compressed canopy (Figure 2A,B), as described by Gou [11]:

$$f_t=f_h\left(1-w\right)+f_{c}w.$$

(1)

Fractions f_h and f_c are calculated by:

$$f_h=1-e^{-k_{cf}LAI},$$

(2)

$$f_c=\left(1-e^{-k_{cf}LA{I}_c}\right)\frac{R}{R+P}.$$

(3)

The leaf area index (LAI) of a compressed canopy LAI_c is defined by:

$$LA{I}_c=LAI\frac{R+P}{R}.$$

(4)

The weighting factor w in Equation (1), which represents the relative contribution of the homogeneous and compressed parts, is defined by:

$$w=\frac{SP-SR}{1-e^{-k_{cf}LA{I}_c}}.$$

(5)

In this equation, SP is the fraction of radiation reaching the path between strips soil surface, and SR the fraction of radiation reaching the soil under the strips:

$$SR= IR{e}^{-k_{cf}LA{I}_c}+\left(1-IR\right)e^{-k_{cf}LAI},$$

(6)

$$SP= IP+\left(1-IP\right)e^{-k_{cf}LAI}.$$

(7)

IP is the fraction which represents the spatial integration of radiation reaching the path between strips, and IR the fraction of representing the spatial integration of the radiation reaching the strip [15]. IP and IR are dependent on the crop plantation form (P and R, respectively) and plant height (H). IP and IR are view factors of the path and plant strip, respectively, and can be calculated according to [26]:

$$IP=\frac{\sqrt{H²+P²}-H}{P},$$

(8)

$$IR=\frac{\sqrt{H²+R²}-H}{R}.$$

(9)

To calculate IP and IR in an intercropping system, the plant height of the taller plant species H_t is used to substitute H in Equations (8) and (9).

To calculate the radiation interception fraction of intercropping plants with different heights, the intercropping canopy is divided into two layers limited by the shorter (H_s) and the taller (H_t) plant height (Figure 2C). The upper layer is limited at the top by H_t and at the bottom by H_s. The lower layer is limited at the top by H_s, and at the bottom by the soil surface. The lower layer is composed by both crop strips alternated. The LAI of the taller canopy is separated into an upper part (Equation (10)), referring to the upper layer of the canopy, and a lower part (Equation (11)) belonging to the lower layer of the canopy.

$$LA{I}_{t,u}=\left(1-\frac{H_s}{H_t}\right)LA{I}_t,$$

(10)

$$LA{I}_{t,l}=\left(\frac{H_s}{H_t}\right)LA{I}_t.$$

(11)

In these equations, LAI_t,u is the portion of the LAI of the taller plant in the upper layer, and LAI_t,l is the portion of the LAI in the lower layer. The compressed form of Equations (10) and (11) can be obtained similarly, replacing the homogeneous plant LAI by the respective LAI for the compressed canopy (LAI_c).

The fraction of the radiation intercepted by the taller crop in the intercropping (Figure 3B) is the sum of radiation intercepted by the taller plant in the upper layer (f_t,u) and in the lower layer (f_t,l) of the canopy:

$$f_t=f_{t,u}+f_{t,l}.$$

(12)

The fraction of the radiation intercepted by the taller plant in the upper layer f_t,u is calculated using Equations (1)–(9). The fraction of the radiation intercepted by the taller plant in the lower layer of the canopy is calculated considering plants to have a compressed arrangement in the field. The lower layer of the canopy consists of the lower part of the taller crops sharing space with the shorter plants. The lower part of the taller plants receives the fraction of the radiation transmitted through the upper part of the taller plant strip (SR_t,u):

$$f_{t,l}=S{R}_{t,u}\left(1-e^{-k_{t}LA{I}_{t,l,c}}\right)\frac{R_t}{R_t+P_t},$$

(13)

$$S{R}_{t,u}= IR{e}^{-k_{t}LA{I}_{t,u,c}}+\left(1-IR\right)e^{-k_{t}LA{I}_{t,u}}.$$

(14)

The fraction of the radiation intercepted by the shorter plants is calculated according to:

$$f_s=S{P}_{t,u}\left(1-e^{-k_{s}LA{I}_{s,c}}\right)\frac{R_s}{R_s+P_s},$$

(15)

where SP_t,u is the fraction the directly incident radiation transmitted through the taller plant strips to the soil path:

$$S{P}_{t,u}= IP+\left(1-IP\right)e^{-k_{t}LA{I}_{t,u}}.$$

(16)

In the case of intercropping species with equal heights (Figure 3C), the fraction of the radiation intercepted by each plant strip is calculated according to:

$$f_{cr}=\left(1-e^{-k_{cr}LA{I}_{cr,c}}\right)\frac{R_{cr}}{R_{cr}+P_{cr}}.$$

(17)

For all plant height combinations in the simulations, the fraction of radiation intercepted by the soil (f_g) is calculated as:

$$f_g=1-f_{t1}-f_{t2}.$$

(18)

2.1.3. Atmospheric Vapor Fluxes

Potential plant transpiration T_p (cm) and soil evaporation E_p (cm) are predicted in the SWAP model using the Penman–Monteith equation [27]. Plant species growing together affect each other’s T_p and E_p. The fraction f_t,i is plant-dependent; therefore, plant transpiration is calculated for each crop (T_p,i) in the simulation unit according to the Penman–Monteith equation with the effect of the presence of the other plant species in the field accounted for by f_t,i in T_p, and the fraction of the radiation intercepted by the soil (f_g) controls the radiation and the aerodynamic terms [28]:

$$T_{p,i}=\frac{\left(1-W_c\right)\left(f_{t,i}\frac{{\Delta }_v}{{\lambda }_w}\left(R_n-G\right)+\frac{p_1{\rho }_a{C}_a}{{\lambda }_w}\left(\frac{e_{sat}-e_a}{r_{a,c,i}}\right)\right)}{{\Delta }_v+{\gamma }_a\left(1+\frac{r_{s,min,i}}{r_{a,c,i }LA{I}_{eff}}\right)},$$

(19)

$$E_p=\frac{f_g\frac{{\Delta }_v}{{\lambda }_w}\left(R_n-G\right)+\frac{p_1{\rho }_a{C}_a}{{\lambda }_w}\left(\frac{e_{sat}-e_a}{r_{a,s}}\right)}{{\Delta }_v+{\gamma }_a\left(1+\frac{r_{soil}}{r_{a,s}}\right)},$$

(20)

The aerodynamic terms of the Penman–Monteith equation (r_a,c,i and r_a,s) are calculated according to:

$$r_{a,c,i}=\frac{r_{a,c,0}}{f_{t,i}},$$

(21)

$$r_{a,g}=\frac{r_{a,s,0}}{f_g}.$$

(22)

In the SWAP model, the actual transpiration T_a,i (cm) is calculated as a function of the soil root distribution and T_p,i, and the combined effect of water, oxygen and salt stress are accounted [24]. In SWAP 2×1D, the assumption of the roots overlapping in the intercropping [29] is not considered and the root distribution model is the SWAP default. The water interception P_in,i (cm) by the canopy of each plant in the intercropping unit is obtained by the method of Von Hoyningen-Hüne and Braden [24]. To distribute P_in,i between plant species canopy, the f_t,i factor of each crop type is added to the P_in,i calculation according to:

$$P_{in,i}=aLA{I}_i\left(1-\frac{1}{1+\frac{f_{t,i}{P}_{gross}}{aLA{I}_i}}\right).$$

(23)

2.1.4. Lateral Soil Water Flux

To calculate the lateral exchange of water between both intercropping strips (q_l,_d, cm d⁻¹) in SWAP 2×1D, the Darcian flow based on the pressure head gradient between two equal-level compartments and their average soil hydraulic conductivity on a daily basis is considered:

$$q_{l,d}=-K_{min}\left(\frac{h_1-h_2}{D_n}\right).$$

(24)

The lateral water balance component is calculated at the end of each simulated 24 h period and then added to the water balance of both SWAP 2×1D modules. SWAP uses a variable but smaller time step (between 1.0∙10⁻⁶ and 0.8 d) to calculate the soil water fluxes and solve the Richards equation, and the addition of a daily water flux at a relatively large time step may cause numerical instability or lead to non-convergence of the numerical scheme in the subsequent time step. To avoid this, a limiting lateral water flux q_lim,d is defined as:

$$q_{lim,d}=\frac{D_n\Delta \theta }{2}.$$

(25)

The actual lateral flux per compartment (q_d) applied in SWAP 2×1D is calculated based on q_lim,d and q_l,d according to:

$$q_d=\frac{q_{lim,d} q_{l,d}}{q_{lim,d}+ q_{l,d}}.$$

(26)

The lateral flux q_d is adjusted to the plant-specific below ground coverage B_c,i (Equation (27)). B_c,i is the proportion of land coverage by the plant species intercropped and its value represents relatively the plant root dominance in the soil profile. The coverage coefficient B_c,i is an input value specified by the user, and the B_c,i values of both plants intercropped must sum up to 1. The relative plant root dominance in soil influences the exchange of water between SWAP 2×1D modules:

$$q_{d,i}=\frac{q_d}{B_{c,i}}.$$

(27)

The total horizontal soil water flux Q_L,i (cm d⁻¹) is calculated as the daily sum of compartment lateral fluxes q_d,i:

$$Q_L={\displaystyle \sum }_{i=1}^n{q}_{d,i}.$$

(28)

2.1.5. Water Balance

The water balance is initially independently calculated for both strips as:

$$\Delta S_1=P-T_{a,1}-E_{a,1}-Q_{b1}-P_{in,1}-R_1+Q_{L,1},$$

(29)

$$\Delta S_2=P-T_{a,2}-E_{a,2}-Q_{b2}-P_{in,2}-R_2-Q_{L,2},$$

(30)

where indices 1 and 2 refer to strips with crops 1 and 2.

The water balance is subsequently calculated for the intercropping simulation unit in SWAP 2×1D from both modules 1 and 2. The intercropping water balance is equal to ∆S_int (Equation (31)), and the sum of ∆S₁ and ∆S₂ multiplied by the respective values of B_c for strip 1 (B_c,1) and strip 2 (B_c,2):

$$\Delta S_{int}= B_{c,1}\Delta S_1+B_{c,2}\Delta S_2.$$

(31)

2.2. Model Evaluation

To evaluate the performance of SWAP 2×1D, a scenario with maize (Zea mays L.) and soybean (Glycine max L.) was simulated using both the single crop SWAP 1D and the strip intercropping SWAP 2×1D. Some details about the simulated scenarios are presented in

Table 2. Overview of the simulated scenarios.

Scenario	Model	Description	Objective
Maize–soybean (MS) intercropping	SWAP 2×1D	Strip intercropping MS (Figure 4A). Soybean is sown approximately 50 days after maize.	Application of SWAP 2×1D to a strip intercropping system
Maize monocropping (M)	SWAP 1D	Monocropping cultivation of maize (Figure 4B)	To compare the water balance outcomes of monocropping to the intercropping scenario.
soybean monocropping (B)	SWAP 1D	Monocropping cultivation of soybean (Figure 4C)	To compare the water balance outcomes of monocropping to the intercropping scenario.

and Figure 4.

For scenarios MS and M, the maize crop cycle started on January 1 (day 1) and ended on April 22, 2015 (day 112), and for MS and S, soybean was simulated to grow from February 20 (day 51) to July 12, 2015 (day 193). To include the entire maize and soybean cycles in the simulations, the simulated period for all scenarios was from January 1 to July 12 (193 days). The distance between two strips of maize (D_s,maize) was 2.0 m in the intercropping scenarios maize-soybean (MS). D_s, strip (R) and path (P) widths are input values for the radiation distribution module in SWAP 2×1D. The strip width R and soil strip width P were both equal to the average 1.0 m and were kept constant during the cropping seasons. R and P were obtained by multiplying plant coverage (B_c) by the distance between rows (D_s). The number of plant rows is absent in the simulations of SWAP 2×1D or SWAP 1D, but it is accounted for by the LAI. Default plant development parameters were used in the WOFOST crop file (

Table 3. Main plant parameters used in the simulation scenarios for maize and soybean (MS), maize (M) and soybean (S).

Description	Parameter	Maize	Soybean	Unit
Plant maximum height	Hmax	200	80	cm
Reflection coefficient, Albedo	Cref	0.20	0.23	-
Minimum canopy resistance	RSC	131	70	s m−1
Temperature sum from emergence to anthesis	Tsum,ea	1000	-	°C
Temperature sum from anthesis to maturity	Tsum,am	1150	-	°C
Maximum CO2 assimilation rate	Amax,d	35	37	kg ha−1 d−1
Maximum relative increase in LAI	RLAI	0.012	0.010	m2 m2
Extinction coefficient for diffuse visible light	kdif	0.60	0.80	-
Extinction coefficient for direct visible light	kdir	0.75	0.80	-
Light use efficiency	eff	0.45	0.40	kg CO2 J−1
Assimilates conversion efficiency into leaves	Cvl	0.68	0.72	kg kg−1
Assimilates conversion efficiency into storage organs	Cvo	0.67	0.68	kg kg−1
Assimilates conversion efficiency into roots	Cvr	0.29	0.72	kg kg−1
Assimilates conversion efficiency into stems	Cvs	0.66	0.69	kg kg−1
Relative increase in respiration rate with temperature	Rit	2.00	2.00	kg CH2O kg−1 d−1
Relative maintenance respiration rate of leaves	Rml	0.03	0.03	kg CH2O kg−1 d−1
Relative maintenance respiration rate of storage organs	Rmo	0.01	0.017	kg CH2O kg−1 d−1
Relative maintenance respiration rate of roots	Rmr	0.015	0.010	kg CH2O kg−1 d−1
Relative maintenance respiration rate of stems	Rms	0.015	0.015	kg CH2O kg−1 d−1
Maximum relative death rate of leaves due to water stress	Pdl	0.03	0.03	d−1
Critical pressure heads according to Feddes [33]	h1	−10.0	−15.0	cm
	h2u	−25.0	−30.0	cm
	h2l	−25.0	−30.0	cm
	h3h	−400.0	−750.0	cm
	h3l	−500.0	−2000.0	cm
	h4	−10000.0	−8000.0	cm
Interception coefficient, corresponding to maximum interception amount	a	0.25	0.25	cm
Maximum daily increase in rooting depth	Rrd,i	2.20	1.2	cm d−1
Maximum root depth	Rd,m	100	100	cm
Below ground plant coverage	Bc	0.5	0.5	-

). Plant H, R_d, and specific leaf area LA were set based on experimental data reported in the scientific literature [29,30,31]. For maximum rooting depths R_d,_m, 100 cm was used both for maize and soybean. The soil was of sandy-clay texture as described in de Jong van Lier [32], with hydraulic properties in

Table 4. Main soil parameters used in the simulation scenarios.

Description	Parameter	Value	Unit
Soil resistance of wet soil	rsoil	230	s m−1
Soil evaporation coefficient of Black	Cbk	0.45	cm1/2
Residual water content	Өr	0.13	cm3 cm−3
Saturated water content	Өs	0.37	cm3 cm−3
Parameter alfa of the water retention curve	α	0.04	cm−1
Parameter n of the water retention curve	n	1.59	-
Exponent in the hydraulic conductivity function	l	1.2	-
Hydraulic conductivity at saturated conditions	Ksat	26	cm d−1

. The soil profile depth is 112 cm, and the profile is subdivided into four sub-layers: 0–12 cm, 12–32 cm, 32–62 cm, and 62–112 cm. Free drainage of the soil profile was used as the bottom boundary condition.

Weather data were taken from the weather station of the University of São Paulo in Piracicaba, Brazil (22°42′ S, 47°38′ W). The total precipitation during the cropping season was 289 mm. No irrigation was applied in the simulation scenarios.

2.3. Sensitivity Analysis

To evaluate the sensitivity of output variables (V_f) to soil and plant parameters, a sensitivity analysis was performed for the intercropping of maize and soybean (

Table 1. List of symbols with description, units, and first equation in which they occur.

Symbol	Description	Equation	Unit
Radiation
ft	Fraction of radiation intercepted by a strip-canopy	(1)	-
fh	Fraction of radiation intercepted by a homogeneous type canopy	(2)	-
fc	Fraction of radiation intercepted by a compressed type canopy	(3)
LAI	Leaf area index over the whole intercropped area	-	m2 m−2
LAIc	Leaf area index of a compressed canopy	(4)	m2 m−2
MLAIi	Maximum leaf area index over the whole intercropped area		m2 m−2
MLAIm	Maximum leaf area index over the whole monocropping area		m2 m−2
w	Factor expressing the relative homogeneous and compressed contribution to the canopy radiation interception model	(5)	-
SR	Fraction of radiation reaching the soil under the strip	(6)	-
SP	Fraction of radiation reaching the soil path	(7)	-
IP	View factor over the soil path	(8)	-
IR	View factor over the strip	(9)	-
LAIt,u	Leaf area index of the taller plants in the upper canopy layer	(10)	m2 m−2
LAIt,l	Leaf area index of the taller plants in the lower canopy layer	(11)	m2 m−2
LAIs,c	Leaf area index compressed of the shorter plant
ft,l	Fraction of radiation interception by the taller plant in the lower layer	(13)	-
ft,u	Fraction of radiation interception by the taller plant in the upper layer	(1)–(9)	-
SRt,u	Fraction of radiation transmitted through the upper part of the taller plants	(14)	-
fs	Fraction of radiation intercepted by the shorter plants	(15)	-
SPt,u	Fraction of radiation transmitted between and through the taller plant to the path	(16)	-
fcr	Fraction of radiation intercepted by strip plants with the same height	(17)	-
fg	Fraction of radiation intercepted by the soil	(18)	-
ft1	Total fraction of the radiation interception by one plant species in the intercropping	(18)	-
ft2	Total fraction of the radiation interception by a second plant species in the intercropping	(18)	-
kcf	Light extinction coefficient	-	-
ks	Light extinction coefficient of the shorter plant	-	-
kt	Light extinction coefficient of the shorter plant	-	-
kcr	Light extinction coefficient of a intercropped plant	-	-
R	Strip width	-	cm
P	Path width	-	cm
Ht	Plant height of the taller plant in the intercrop	-	cm
Hs	Plant height of the shorter plant in the intercrop	-	cm
Rt	Strip width of the taller plant in the intercrop	-	cm
Pt	Path width of the taller plant in the intercrop	-	cm
Rs	Strip width of the shorter plant in the intercrop	-	cm
Ps	Path width of the shorter plant in the intercrop	-	cm
Rcr	Strip width of a intercropped plant	-	cm
Pcr	Path width of a intercropped plant	-	cm
Water flux
Tp,i	Potential transpiration of each plant separately	(19)	cm
Ep	Potential evaporation	(20)	cm
Wc	Fraction of the day during which the canopy is wet	(19)	-
∆v	Slope of the saturated vapor pressure-temperature curve	(19), (20)	kPa °C−1
λw	Latent heat of vaporization	(19), (20)	kJ kg−1
Rn	Net radiation at the canopy surface	(19), (20)	J m−2 d−1
G	Soil heat flux density	(19), (20)	J m−2 d−1
p1	Unit conversion value	(19), (20)	s d−1
ρa	Air density	(19), (20)	kg m−3
Ca	Heat capacity of the moist air	(19), (20)	J kg−1 °C−1
esat	Saturation vapor pressure	(19), (20)	kPa
ea	Actual vapor pressure	(19), (20)	kPa
γa	Psychrometric constant	(19), (20)	kPa °C−1
ra,c,i	Aerodynamic resistance of the plant canopy	(21)	s m−1
ra,g	Aerodynamic resistance of the soil	(22)	s m−1
rs,min,i	Minimum stomatal resistance	(19)	s m−1
LAIeff	Effective leaf area index	(19)
rsoil	Soil resistance of a wet soil	(20)	s m−1
Pin,i	Water interception by the plant canopy	(23)	cm
Pgross	Gross precipitation	(23)	cm d−1
a	Empirical coefficient for the precipitation interception	(23)
ql,d	Lateral Darcian water flux between soil compartments	(24)	cm d−1
Kmin	Hydraulic conductivity of the driest same-level soil compartment	(24)	cm d−1
Dn	Distance between the centers of two neighboring strips of different plant species	(24), (25)	cm
h1, h2	Soil pressure head	(24)	cm
qlim,d	Limiting lateral flux per compartment	(25)	cm d−1
∆θ	Difference in soil water content of two same-level soil compartments	(25)
qd	Lateral flux per soil compartment limited by qlim	(26)	cm d−1
Bc,i	Plant-ground coverage fraction	(27), (31)	cm2 cm−2
QL	Total lateral flux of the soil profile	(27)	cm d−1
Water balance
Pa	Precipitation	(29), (30)	cm
Ta	Plant transpiration	(29), (30)	cm
Ea	Soil evaporation	(29), (30)	cm
Qb	Bottom flux	(29), (30)	cm
∆S1	Variation of soil water storage in soil-plant module 1	(29), (31)	cm
∆S2	Variation of soil water storage in soil-plant module 2	(30), (31)	cm
∆Sin	Variation of soil water storage for the intercropping	(30), (31)	cm
Sensitivity analysis
η	Relative partial sensitivity index	(32)	-
Vf	Default value of a generic output variable	(32)
Vs	New value of the output variable after a selected parameter was changed	(32)
Vr	Reference value of the output variables evaluated in the sensitivity analysis	(32)
p	General representation for parameter	(32)
∆p	Change in the parameter values in the sensitivity analysis	(32)
Cref	Canopy reflection coefficient		-
RLAI	Maximum relative increase in LAI		m2 m−2 d−1
kdif	Extinction coefficient for diffuse light		-
kdir	Extinction coefficient for direct light		-
Өr	Saturated soil water content		cm3 cm−3
Өs	Residual soil water content		cm3 cm−3
α	Soil water retention curve parameter		cm−1
n	Soil water retention curve parameter		-
Ks	Saturated hydraulic conductivity		cm d−1
l	Soil hydraulic conductivity exponent		-

). The evaluated output variables are: radiation interception fraction by the intercrop (f_t), actual crop transpiration (T_a), actual soil evaporation (E_a), lateral soil water flux (Q_L) and flux below the root zone (Q_b). Each of these output variables were evaluated using the relative partial sensitivity index η:

$$\eta =\frac{\left(V_f-V_s\right)}{V_r}\frac{p}{\Delta p},$$

(32)

in which V_f is the evaluated output variable, V_s is the new value of the output variable after a selected parameter was changed, and (V_f − V_s) the difference of an output variable due to the change of parameter p. The reference V_r is the average of T_a, E_a, Q_L, and Q_b output values for water fluxes, and the average of the output values of f_t for the radiation interception simulated with the default (reference) parameters values. To calculate η, the selected parameter was changed by 1% (Δp/p = 0.01) while the others were maintained at their default value. The reference value V_r was 9.5 cm for the soil fluxes variables T_a, E_a, Q_L, and Q_b, and 0.38 for the radiation interception f_t. A variable with |η| ≥ 0.5 was considered to have high sensitivity to the parameter evaluated.

3. Results and Discussion

3.1. LAI and Radiation Interception

Intercropped maize LAI simulated by SWAP 2×1D was lower than in maize monocropping (Figure 5A) for most of the simulation period. This difference in LAI can be explained since in the intercropping maize–soybean plants of soybean were occupying half of the cultivation area covered by maize plants in the monocropping. Maize LAI was calculated for the whole cultivated area in SWAP 2×1D and for that reason it is lower than the monocropping maize LAI predicted by SWAP 1D. Intercropped soybean plants were simulated to grow together with taller plants of maize and to receive less radiation compared to the monocropping cultivation of soybean, another reason for soybean intercropping LAI to be reduced (Figure 5B). SWAP 2×1D was able to reproduce the reduction in the intercrop LAI compared to the monocropping system, similarly to the behavior of the intercropping LAI observed by Gou [11] and Liu [34].

The fraction of radiation interception (f_t) in maize as a function of time showed similar to the intercropping LAI, showing the low interference of the shorter soybean plants in this scenario. Soybean f_t increased suddenly on day 112 because maize plants were harvested on this date and no longer covered the remaining soybean plants. The maximum f_t by the intercropped maize canopy simulated by SWAP 2×1D was 0.80 at LAI 5.40 m² m⁻² and maximum f_t by soybean canopy was 0.49 at LAI 1.38 m² m⁻² (Figure 5C).

The fraction of radiation intercepted by the soil surface (f_g) suddenly decreased at the end of maize cycle on day 112 (Figure 5D). After maize harvest, more radiation reached the soil, attenuated only by the soybean plants. SWAP 2×1D does not consider the possibility of existing plant shoots and remaining biomass over the soil after harvest, which could result in a smaller value of f_g.

3.2. Transpiration and Evaporation Fluxes

Intercropped maize T_a is higher than maize monocropping T_a for most of the cultivation period (Figure 6A). The maximum T_a for intercropped maize is 1.05 cm d⁻¹ and for maize monocropping is 0.68 cm d⁻¹ at the end of the vegetative cycle (day 55). The effect of shadowing of soybean by maize is observed between days 50 and 100, when soybean (MS) T_a is 0.05 cm d⁻¹ and monocropping soybean (S) T_a is 0.14 cm d⁻¹. The maximum T_a for intercropped soybean is 0.30 cm d⁻¹ and for monocropping soybean is 0.41 cm d⁻¹ (on day 92, Figure 5B). The sudden reduction in soybean monocropping and intercropping T_a after day 95 is a response to water shortage, however. On average, intercropped soybean T_a is reduced by 0.05 cm d⁻¹ and monocropping soybean T_a by 0.06 cm d⁻¹ due to water shortage from day 95 to 192.

The transpiration fluxes obtained in these hypothetical scenarios were not validated, but the intercropped plant transpiration behavior (T_a) simulated by SWAP 2×1D in comparison to their monocropping T_a values was in agreement with the behavior of the intercropping wheat-maize and wheat-sunflower compared to their monocropping systems as obtained by Miao [31].

Intercropped maize (MS) soil E_a predicted by SWAP 2×1D was correspondent to maize monocropping (M) soil E_a in most of the season except a few periods (Figure 6C). From day 50 to 80, differences in MS and M E_a are probably an effect of distinct plant transpiration in both scenarios. Between days 150 to 200, the soil was uncovered in scenario M, while in MS the soil evaporation was affected by soybean strips and their partial interception of the incoming radiation, causing the difference in evaporation between M and MS. The intercrop effect of maize plants on soybean E_a are simulated by SWAP 2×1D (Figure 6B,D). From day 41 to 80, in the intercropping period, average E_a of the soil under intercrop soybean is 0.04 cm d⁻¹ while E_a in monocropping soybean is 0.20 cm d⁻¹ (Figure 6D).

3.3. Lateral Soil Water Flux

The higher water demand by maize plants made SWAP 2×1D predict a lateral water flow (average 0.15 cm d⁻¹) from the soybean towards the maize strips (Figure 7A) during most of the time until day 112. In the first 20 days of maize development, water was supplied by the frequent rainfall events and no significant lateral flow developed. Intercropped maize water demand increased, making the lateral water flow to increase from day 20 to 48. As Q_L is governed by soil pressure head gradients as well as by K_min, peak events of Q_L (e.g., on days 47, 48, 56 and 151) were caused by events of high precipitation causing an instantaneous increase of K_min (Figure 7A,B). After the end of the maize cycle (day 112), water flows laterally at an average 0.03 cm d⁻¹ towards the soybean soil region until the end of the soybean crop season.

The significance of the lateral water flow for soil water availability in the intercropping (MS) is observed comparing the soil water content profiles of maize and soybean predicted by SWAP 2×1D with and without the horizontal soil water exchange (Figure 8). Within the first 48 days, maize strips are cultivated in between paths of uncovered soil (Figure 8A), and the soil water content (θ) in the 55–65 cm depth region in maize soil profile is greater than in maize (MS) simulation with no lateral flux. The θ of the soil path (strip of uncovered soil) in the soil depth 0–80 cm for the simulations with lateral water flow is lower than the θ of the soil profile in the simulations without the lateral water flow. On day 71, the effect of the lateral soil water flow on θ of the intercropping soil profile is evident (Figure 8B). Soybean soil θ distribution tends to be lower in the simulation with lateral flux since water was flowing preferentially from soybean to maize until day 71.

3.4. Water Balance

The water balance components of the intercropping MS and monocropping maize (M) and soybean (S) are shown in

Table 5. Water balance components for intercropping MS and monocropping maize and soybean: ∆S, water storage change; P, precipitation; Q_b, bottom flux; P_i, canopy water interception; T_a, actual transpiration; E_a, actual evaporation; ET_a and ET_P, actual and potential evapotranspiration, respectively; T_p, potential transpiration; E_p, potential evaporation.

Water Balance Components	SWAP 1D			SWAP 2×1D
(cm)	Maize	Soybean	MS 1D Average	MS
Actual
∆S	4.1	1.1	2.6	1.1
Pa	47.7	47.7	47.7	47.7
Pi	−2.7	−1.0	−1.9	−1.6
ETa	−40.8	−37.2	−39.0	−41.6
Ta	−24.6	−17.7	−21.1	−25.4
Ea	−16.2	−19.5	−17.9	−16.1
Qb	−0.1	−8.3	−4.2	−3.5
Potential
ETp	74.3	75.0	74.7	73.2
Tp	24.6	24.2	24.4	31.4
Ep	49.7	50.8	50.3	41.8

. The canopy water interception (P_i), soil evaporation (E_a), and bottom flux (Q) for the intercropping simulated with SWAP 2×1D were similar to the average of the monocropping water balance components simulated with SWAP 1D. MS intercropping potential and actual transpiration rates (T_p and T_a) predicted by SWAP 2×1D are higher than the weighted averages T_a and T_p of the monocropping maize and soybean obtained with SWAP 1D. Consequently, other components of the water balance, the water storage (S) and Q_b for the intercropping MS, are reduced in comparison to average of the monocropping maize and soybean. The actual and potential evapotranspiration (ET_a and ET_p) were similar for the scenarios M, S, simulated with SWAP 1D, and MS with SWAP 2×1D.

The lateral water flux (Q_L) simulated with SWAP 2×1D for the whole period of MS intercropping is 5.0 cm, from the soybean to the maize soil strip. Q_L is significantly low compared to the other components of the water balance ET_a, T_a and E_a. However, the low exchange of water laterally observed is the resultant sum of the daily positive (soybean to maize) and negative (maize to soybean) lateral fluxes. Positive lateral fluxes Q_L in the soil under maize occurred frequently in the first half of the whole simulation period (Figure 7). With the end of maize cycle, the lateral water flux was frequently in the soybean direction. The maximum Q_L was 4.1 cm d⁻¹, and the water flow from maize to soybean strips. The wet period during maize-strips and MS intercropping is also predominant for the resultant low lateral flux. A scenario of intercropping in which water is not available sufficiently and plants are sown at the same time can result in distinct lateral flow behavior and resultant Q_L.

3.5. Sensitivity Analysis

Intercropped maize potential transpiration (T_p) is sensitive to maize radiation use efficiency (RUE), albedo (C_ref), maximum increase of LAI (RLAI), and coefficients of light extinction (k_dif and k_dir). Since T_p is directly related to the radiation interception by the plant, an increase in C_ref results an increase in the radiation reflection and decrease in the radiation absorption and, consequently, a reduction in T_p. An increase in RUE, R_LAI, k_dif, and k_dir increased T_p due to an improvement in the conversion of the incoming radiation.

Maize plant parameters albedo (C_ref) and maximum increase of LAI (R_LAI) show significant influence on maize actual transpiration (T_a), lateral flux (Q_L), and bottom flux (Q_b). Accordingly, changes in maize C_ref and R_LAI affect intercropped soybean Q_L and Q_b significantly. Intercropped maize Q_L and Q_b are sensitive to C_ref and R_LAI, as a consequence of T_a behavior in the sensitivity analysis.

Maize T_a is affected by the saturated soil water content (Ө_s), shape parameters of the soil water retention curve a and n, and the saturated hydraulic conductivity (K_s), which are the parameters that change effectively the soil water availability for plants in SWAP. The lateral flux Q_L and bottom flux Q_b were sensitive to most of the evaluated soil parameters.

The variables presented in

Table 6. Relative partial sensitivity index η of soil and atmospheric fluxes for plant and soil parameters predicted by the SWAP 2×1D model.

Parameter	η
	Maize							Soybean
	∑ft	Tp	Ta	Ep	Ea	QL	Qb	∑ft	Tp	Ta	Ep	Ea	QL	Qb
Maize
RUE	0.17	0.63	−0.11	0.00	−0.21	−0.21	−0.42	0.00	0.00	0.00	−0.42	−0.11	0.21	0.00
Cref	0.01	−1.16	1.05	0.00	0.00	1.05	2.21	−0.01	0.00	0.00	0.00	0.00	−1.05	−1.16
RLAI	0.27	1.26	−2.00	−0.01	−0.42	−2.31	−4.62	−0.01	0.00	0.11	−0.74	−0.21	2.31	2.00
kdif	0.19	0.53	0.00	0.00	−0.42	−0.42	−0.95	−0.07	−0.21	0.11	−0.42	−0.11	0.42	0.21
kdir	0.40	1.26	0.42	−0.01	−0.53	−0.11	−0.11	−0.06	−0.21	0.11	−0.84	−0.21	0.11	−0.32
Soil
Өr	0.00	0.00	0.32	0.00	0.11	0.74	1.47	−0.18	−0.42	0.63	0.21	−0.11	−0.74	−1.58
Өs	0.01	0.00	5.78	0.00	1.26	7.04	14.18	0.09	−0.32	0.53	0.11	1.05	−7.04	−7.25
a	−0.03	−0.11	−1.16	0.00	0.00	−0.74	−1.47	−0.12	−0.32	0.42	0.21	0.00	0.74	0.53
n	0.00	0.00	−1.37	0.00	−0.42	−0.53	−1.16	−0.16	−0.32	0.63	0.21	−0.21	0.53	−0.21
Ks	0.00	0.00	1.47	0.00	0.11	1.79	3.68	0.01	0.00	0.00	−0.11	0.00	−1.79	−2.31
l	0.00	−0.11	−0.21	0.00	−0.11	−0.21	−0.42	0.02	0.00	0.00	0.00	0.00	0.21	0.32

Note: Output variable: ∑f_t, the sum of daily plant interception fractions; T_p and T_a, potential and actual transpiration rates; E_p and E_a, potential and actual evaporation rates; Q_L-in, lateral influx; Q_b, bottom flux. Parameter: RUE, radiation use efficiency; C_ref, reflection coefficient (Albedo); R_LAI, maximum relative increase in LAI; k_dif and k_dir, extinction coefficients for diffuse and direct light, respectively; Ө_r and Ө_s, the residual and saturated soil moistures, a and n, retention curve parameters, l, soil hydraulic conductivity parameter, and K_s, the saturated hydraulic conductivity.

had no significant sensitivity to soybean parameters in the evaluation, and their results were, therefore, omitted. The highest sensitivities (η > 2) of maize and soybean simulations were obtained for Ө_s and R_LAI.

Most of the plant and soil parameters included in Table 6 are important for the simulation of intercropping MS. These parameters are recommended to be adjusted in a calibration and validation processes of the model SWAP 2×1D. Nevertheless, the sensitivity of the output variables available in Table 6 refers to SWAP 2×1D simulated with the default parameters and meteorological conditions in the current hypothetical scenario. The sensitivity of the output variables to the selected parameters may be different in a distinct scenario.

4. Conclusions

Radiation and soil water sharing subroutines were implemented in the model SWAP/WOFOST 1D for the simulation of strip intercropping cultivation. We presented theory and concepts used in the newly developed model SWAP 2×1D and simulated hypothetical scenarios to discuss the significance of simulating the lateral soil water flux in the intercropping in addition to the radiation sharing between plants. SWAP 2×1D was capable of simulating above- (light) and belowground (water) interactions in strip intercropping.

The lateral water flow between soil strips in the intercropping scenario was relatively low compared to other components of the water balance on a cumulative basis. Daily values, however, were of the same order of magnitude as transpiration and evaporation rates. The effect of the lateral water exchange on the soil water storage of the intercropping was shown to be significant.