A method for improving soil and fixing carbon by coordinating straw returning to field and tillage sequence

Abstract

The present application relates to the technical field of soil improvement, and discloses a method for soil improvement and carbon sequestration by synergistically regulating straw returning and tillage sequences. Corn stalks are classified into fine particle size modes and coarse particle size modes with non-overlapping particle size intervals and are compounded as conditioning materials, and marked conditioning materials and unmarked conditioning materials with different carbon-13 atomic fractions but the same composition are prepared; shallow strips and deep strips with different soil layer depths are opened in the same field and filled with the conditioning materials, and shallow marked sub-regions, deep marked sub-regions and reference sub-regions are divided; multiple carbon pools are separated by stratified sampling before deep tillage, and the mineral combination recovery rate and the closed storage protection recovery rate are calculated based on the background of the same carbon pool in the reference sub-region; the phase coordinate target value is calibrated by a response surface model, and the deep tillage is arranged in the 1st basic tillage operation after both the phase coordinate target values are reached. The present application improves the net carbon sequestration and the stability of soil aggregates under the condition of the same number of deep tillage operations.

Description

Technical Field

[0001] This invention relates to the field of soil improvement technology, and in particular to a method for soil improvement and carbon sequestration through synergistic regulation of straw return to the field and tillage sequence. Background Technology

[0002] Returning straw to the field is a major way to improve soil organic matter in farmland. Current technologies typically involve rotary tilling and burying crushed straw in the topsoil, or using deep loosening and deep plowing to push straw into deeper soil layers, using the amount of straw returned, the depth of return, the tillage depth, and the frequency of tillage as control variables. Existing studies have used carbon-13 tracers to compare the distribution of straw carbon between particulate and mineral-bound states under different burial depths. Other studies have laid straw as a layer and observed its impact on underlying organic carbon, indicating that the depth of return alters the fate of straw carbon. However, several problems remain with the control methods used in production. Typically, only one return depth can be set within a single field. The depth effect is intertwined with the input amount effect, soil aeration, and differences in hydrothermal environment, making it impossible to trace the fate of carbon input at a particular depth in isolation. Tillage control generally uses intensity and frequency as variables. Even with one deep tillage operation within a crop rotation cycle, the year in which this operation takes place is often treated as a matter of agricultural timing, lacking a quantitative description of its impact on carbon sequestration. The timing of deep tillage is mainly determined by experience, lacking means to assess the progress of straw-derived carbon conversion before operations. If disturbance occurs during the peak period of carbon formation in the closed aggregate state, the physically protected straw-derived carbon will be re-exposed and mineralized more rapidly. Carbon sequestration effects are mostly characterized by changes in organic carbon content at a fixed soil depth, which neither accounts for the activation loss of primary organic carbon due to the input of exogenous materials nor takes into account the impact of mechanical operations and nitrous oxide emissions, thus limiting the reliability of the conclusions. Summary of the Invention

[0003] This invention aims to address the problems in existing straw return and tillage technologies, such as the difficulty in tracing the destination of carbon input at different straw return depths within the same field, the lack of measurable criteria for determining the timing of deep tillage operations, and the incomplete calculation of carbon sequestration effects. It provides a soil improvement and carbon sequestration method that synergistically regulates the straw return and tillage sequence.

[0004] The technical solution of this invention is implemented as follows: A method for soil improvement and carbon sequestration that synergistically regulates straw return to the field and tillage sequence includes the following steps: Corn stalks were graded into fine and coarse particle size modes with non-overlapping particle size ranges and compounded into conditioning materials according to a set mass ratio. Labeled and unlabeled conditioning materials with different carbon 13 atom fractions but the same particle size mode composition and chemical composition were prepared. Within the same field, shallow and deep strips with different soil depths are cut according to the set number of strips and filled with conditioning material. At the same time, the field is divided into three sub-zones. The shallow marked sub-zone is filled with marked conditioning material only in the shallow strips, the deep marked sub-zone is filled with marked conditioning material only in the deep strips, and the reference sub-zone is filled with unmarked conditioning material. The application parameters for the three sub-zones are the same as those for field operations. Before deep tillage, soil samples were collected from each sub-region in layers to separate multiple non-overlapping carbon pools, including closed-state granular organic carbon pools and fine-grained heavy component organic carbon pools. Using the carbon 13 atomic fraction of the carbon pool of the same name in the reference sub-region as the baseline, the labeled source carbon reserves of each carbon pool were calculated based on the carbon 13 atomic fraction. The ratio of the labeled source carbon reserves of the fine-grained heavy component organic carbon pool and the closed-state granular organic carbon pool to the labeled carbon input was used as the mineral binding recovery rate and the closed-state protection recovery rate, respectively. The target value of the phase coordinate, which is composed of the mineral binding recovery rate and the closed storage protection recovery rate, was obtained by calibration using the response surface model. Deep tillage was scheduled to be carried out in the first basic tillage operation after the mineral binding recovery rate and the closed storage protection recovery rate reached the target value of the phase coordinate. Rotary tillage was used for the remaining basic tillage operations in the rotation cycle.

[0005] Furthermore, the particle size of the fine-particle-size mode is 0.5 mm to 2 mm, the length of the coarse-particle-size mode is 20 mm to 50 mm, and a particle size gap of 2 mm to 20 mm is set between the two modes; the mass ratio of the coarse-particle-size mode to the fine-particle-size mode is 3:1, 1:1, or 1:3; the carbon-13 atomic fraction of the labeled conditioning material is 0.0165 to 0.0200, and the carbon-13 atomic fraction of the unlabeled conditioning material is 0.01080 to 0.01092; the relative deviations between the labeled and unlabeled conditioning materials in terms of organic carbon mass fraction, carbon-nitrogen ratio, lignin mass fraction, cellulose mass fraction, hemicellulose mass fraction, and water-soluble organic carbon mass fraction under the same particle size mode are all no greater than 5%.

[0006] Furthermore, the labeled and unlabeled processed materials are harvested from two cultivation chambers of the same maize variety that are cultivated in parallel in two cultivation chambers with different carbon 13 atomic fractions of carbon dioxide but the same other cultivation conditions. There are no fewer than three cultivation batches, and the carbon 13 atomic fractions of carbon dioxide supplied by the two cultivation chambers are exchanged between adjacent batches.

[0007] Furthermore, the shallow strips are located in the 10cm to 15cm soil layer, and the deep strips are located in the 25cm to 30cm soil layer. The spacing between adjacent strips is 0.47m, the strip width is 0.09m, the strip thickness is 0.025m, the bulk density is 140kg / m3 to 160kg / m3, and the carbon loading per unit strip length is 0.135kg / m to 0.150kg / m. The carbon input per unit area is 2900kg / hm2 to 3200kg / hm2. The ratio of the number of shallow strips to deep strips in every 4 consecutive strips is 3:1, 2:2, or 1:3.

[0008] Furthermore, the area was divided into double-marked sub-areas and blank sub-areas. In the double-marked sub-areas, both shallow and deep strips were filled with marked conditioning material. In the blank sub-areas, trenches were dug according to the same strip ratio and backfilled with an equal volume of original soil. The sub-areas were randomly arranged, with an isolation zone at least one machine operating width between adjacent sub-areas. The marked conditioning material was laid out along the strips in 0.20m long marked sections, directly contacting the soil, with a minimum distance of 0.50m between adjacent marked sections. Each marked section had a 0.15m long unmarked conditioning material protection section at both ends. Sampling and excavation were then carried out. The length is not less than the sum of the length of the marked segment and twice the maximum lateral migration distance, which is determined by the inert tracer pre-test; the relative deviation between the carbon source reserves of the two-position marked subregion and the sum of the carbon source reserves of the shallow-position marked subregion and the deep-position marked subregion is not greater than 10% as an additivity check; the carbon source ratio of the material is calculated by dividing the difference of the carbon 13 atomic fraction between the reference subregion and the blank subregion by the difference of the carbon 13 atomic fraction between the unlabeled conditioning material and the blank subregion carbon library of the same name as a diagnostic check.

[0009] Furthermore, the carbon pool separation and source carbon storage calculation were performed as follows: Soil samples were sieved through a 2mm sieve. Identifiable residues on the sieve were picked out, rinsed, dried, weighed, and directly measured to obtain the coarse residue carbon pool. Soil samples below the sieve were added to a sodium polytungstate solution with a density of 1.8 g / cm³ and centrifuged at 1000 times gravity acceleration for 15 min. The floating matter was the free particulate organic carbon pool. The sediment was ultrasonically dispersed at an energy density of 300 J / mL and then subjected to repeated density separation. The second floating matter was the closed-state particulate organic carbon pool. The remaining heavy components were separated by 0.02 m... The soil samples were separated into fine-grained heavy organic carbon libraries and coarse-grained heavy organic carbon libraries, with m as the boundary. Inorganic carbon was removed from carbonate-containing soil samples by hydrochloric acid vapor fumigation. The mass recovery rate of the five carbon libraries was greater than 95%, and the carbon recovery rate was 95% to 105%. The carbon 13 atomic fraction of each carbon library was obtained by dividing the difference between the carbon 13 atomic fraction of the labeled conditioner and the carbon 13 atomic fraction of the same carbon library in the reference subregion by the difference between the carbon 13 atomic fraction of the labeled conditioner and the carbon 13 atomic fraction of the same carbon library in the reference subregion. The labeled source carbon ratio of each carbon library was then calculated from the labeled source carbon ratio, the organic carbon content of the carbon library, the soil bulk density, and the soil thickness.

[0010] Furthermore, under no-till conditions, samples were taken once each on days 0, 7, 30, 90, 180, 360, 720, 1080, and 1440 after the addition of conditioning materials. Hydrothermal time was accumulated according to soil depth in the strips: data from 10cm to 15cm soil layer for shallow strips and data from 25cm to 30cm soil layer for deep strips. Hydrothermal time was the cumulative difference between the daily average temperature and 5℃ for each day meeting the requirement of a daily average temperature above 5℃ and soil moisture content not less than 50% of field capacity. A model was established for the transformation of coarse residual carbon pools into free particulate organic carbon pools, and the transformation of free particulate organic carbon pools into closed-system particulate organic carbon pools, fine-grained heavy component organic carbon pools, and coarse-grained heavy component organic carbon pools. The first-order kinetic equations, with rate terms representing the transformation of the carbon pool, the transformation of the closed-state particulate organic carbon pool and the coarse-grained heavy component organic carbon pool to the fine-grained heavy component organic carbon pool, and the losses of each carbon pool in the form of mineralization and leaching, were used as the initial values ​​of the measured values ​​on day 0. A nonlinear mixed-effect model was established and fitted using the natural logarithm of each rate constant as a parameter. Before fitting, the identifiability of the parameters was tested by the rank and condition number of the sensitivity matrix. The rank-deficient parameters were fixed one by one according to the Akaike information content criterion. The normalized trajectories of the shallow and deep strips were synthesized by weighting them according to the proportion of their respective labeled carbon inputs, with the calendar time as a common index. The subsequent crop was sown on the day when the cumulative hydrothermal time in the 0cm to 20cm soil layer reached 240℃·d, and the latest sowing date was set.

[0011] Furthermore, a full crossover experiment was conducted with three factors: the ratio of conditioning materials, the proportion of strips, and the tillage sequence. A total of 27 treatments were implemented, repeated three times, and started in three consecutive years, resulting in a total of 81 plots. Alternatively, a D-optimal design was adopted, which allowed for the estimation of all interaction terms between the two factors and had no fewer than 24 different design points. The three tillage sequences were all deep-tilled once within a four-year crop rotation cycle, but in different years. The soil depth for deep-tillage was 30 cm, and the soil depth for rotary tillage was 15 cm. A blank subregion with the same proportion of strips was set up within each plot.

[0012] Furthermore, prior to the experiment, soil columns ranging from 0cm to 60cm were collected from each plot, and the bulk density, cumulative soil mass, percentage of particles smaller than 0.02mm, organic carbon content in each carbon pool, and percentage of water-stable aggregates larger than 0.25mm were measured in 5cm layers as baseline values. Net carbon sequestration was calculated as follows: the minimum cumulative soil mass from 0cm to 45cm across all plots was used as the baseline soil mass, and organic carbon storage was calculated using the equal-mass soil layer method. The increment of organic carbon storage in the treated plot relative to the baseline value and the increment of organic carbon storage in the blank sub-region relative to the baseline were taken. The difference in the increments of values ​​was converted into carbon dioxide equivalents at a ratio of 44 to 12, and then the carbon dioxide equivalents of the cumulative increase in methane emissions converted at 27, the carbon dioxide equivalents of the cumulative increase in nitrous oxide emissions converted at 273, and the emissions of the increase in diesel consumption for tillage operations converted at 2.68 kg of carbon dioxide per liter were subtracted. Greenhouse gas fluxes were measured every 7 to 10 days throughout the growing season, and were measured more frequently for 5 consecutive days after tillage operations, fertilization, daily precipitation greater than 10 mm, irrigation, and freeze-thaw events. Soil carbon dioxide fluxes were counted separately and used to verify mineralization and activation processes.

[0013] Furthermore, the target values ​​for the phase coordinates are calibrated as follows: Net solid carbon is used as the dependent variable, and the independent variables are the linear and quadratic terms of mineral-bound recovery rate and closed-loop protection recovery rate, their interaction term, the main effects of conditioning material mass ratio and strip number ratio, their interaction term, carbon loading per unit strip length, and the random intercept of the start-up year. A response surface model is fitted, with the tillage sequence not included as a factor in the response surface model. Before fitting, a coordinate identifiability test is performed. The coordinate identifiability test is performed if the variance inflation factors of both mineral-bound recovery rate and closed-loop protection recovery rate are less than 2, the design matrix is ​​full rank and the condition number is less than 30, and the constraint optimum lies within the convex hull formed by the measured values ​​of mineral-bound recovery rate and closed-loop protection recovery rate. The total retention rate is defined as the ratio of the sum of carbon reserves from marked sources in each carbon pool to the marked carbon input, with a 95% confidence limit of not less than [value missing]. Under the constraints of a pre-set threshold and a water-stable aggregate mass ratio higher than the baseline value, the optimal value of the response surface model is determined to obtain the target value of the phase coordinate. An alternative model with the tillage sequence as a factor is then fitted, and the alternative model is compared with the response surface model according to the Akaike information content criterion and cross-validation. When the mineral binding recovery rate and the closed-loop protection recovery rate do not simultaneously reach the target value of the phase coordinate within one crop rotation cycle, rotary tillage is adopted for all basic tillage operations within the same crop rotation cycle. The number of deep strips is selected based on the largest net carbon fixation amount predicted by the response surface model, under the constraint that the fine-grained heavy component organic carbon saturation deficit is not less than 0. The saturation deficit is the difference between the estimated value of the soil fine-grained heavy component organic carbon capacity and the measured fine-grained heavy component organic carbon content, both expressed as the mass of carbon contained in each kilogram of soil. The final output includes the conditioner composition, strip application parameters, target value of the phase coordinate, and tillage sequence.

[0014] This invention offers the following advantages: Both labeled and unlabeled conditioning materials are derived from the same batch of raw materials, and their particle size modal composition and chemical composition are consistent except for the carbon-13 atomic fraction. Furthermore, the application parameters and field operations for each sub-region are identical. Therefore, the destination of carbon input at different return depths can be traced within the same field, separating the effect of return depth itself from the effects of input amount, soil environmental differences, and the activating effects caused by organic material input. By using the same carbon pool as the reference sub-region and other sub-regions to deduct background disturbances, the interferences of background organic carbon, current crop root carbon and rhizosphere deposited carbon, and the redistribution of background carbon caused by cultivation can be simultaneously eliminated. The timing of deep tillage operations has been changed from relying on experience to determine based on measured phase coordinates composed of mineral-bound recovery rate and closed-cell recovery rate. Mechanical disturbance thus occurs during periods when a significant portion of straw-derived carbon has been transferred to the mineral surface and the portion relying on physical agglomeration has receded. This allows for an increase in net carbon sequestration and the proportion of water-stable agglomerates without altering the number of deep tillage operations or the amount of straw input. The ratio of coarse to fine mass in conditioning materials can change the shape of the transformation trajectory, thereby altering the timing of reaching the target phase coordinate. Material composition and tillage sequence can thus be jointly optimized. Net carbon sequestration is calculated using the double-difference method and the equal-mass soil layer method, incorporating emission increments. The calculation results accurately reflect the net change in the soil carbon pool. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the soil profile principle for the application of fixed-depth strips and the replacement coding of strip labels according to the present invention; Figure 2 This is a planar schematic diagram of the sub-region coding layout, marking segments, protection segments, and sampling excavation range of the present invention. Figure 3 This is a schematic diagram illustrating the transformation pathways of source carbon among five carbon pools and the principle of tillage sensitivity zoning in this invention. Figure 4 The graph shows the transformation trajectory of the carbon reserves from the source marked by the present invention as a function of calendar time, wherein (a) is the transformation trajectory of the shallow stripe and (b) is the transformation trajectory of the deep stripe. Figure 5 This is a schematic diagram of the phase coordinate trajectory and the target phase coordinate value of the present invention; Figure 6 This is a schematic diagram of the contour lines of the response surface of net solid carbon content to phase coordinates and the constraint optimal point of the present invention. Detailed Implementation

[0016] The embodiments of the present invention will be further described below with reference to the accompanying drawings. The present invention deals with farmland soil in a two-crop-a-year rotation of wheat and corn, with a four-year rotation cycle. Basic tillage is performed once after the corn harvest each year, for a total of four basic tillage operations within one rotation cycle. Basic tillage operations are divided into two types: deep tillage with a soil depth of 30cm and rotary tillage with a soil depth of 15cm. The coded strips are laid once after the corn harvest in the first year of the rotation cycle, and the subsequent four basic tillage operations fall on the 365th, 730th, 1095th, and 1460th days after laying the strips. The present invention does not change the number of deep tillage operations within one rotation cycle, nor does it change the total amount of straw input. Instead, it changes the year sequence of deep tillage operations from judgment based on agricultural experience to determination based on the measured conversion progress of straw-derived carbon in the soil, and simultaneously determines the corresponding material composition and strip application parameters.

[0017] After entering the soil, straw undergoes several states of existence: fresh straw initially exists as identifiable remnants, then breaks down into free particulate organic matter. A portion of this organic matter is encapsulated by newly formed aggregates, while another portion is transformed into organic carbon bound to fine minerals through microbial assimilation and adsorption to mineral surfaces. These two stabilization pathways respond completely differently to mechanical disturbance: aggregate encapsulation relies on the integrity of the aggregates themselves; deep tillage breaks up these aggregates, exposing the encapsulated organic matter to microorganisms and oxygen, significantly increasing the mineralization rate. Organic carbon bound to fine mineral surfaces, on the other hand, is fixed by chemical bonding and adsorption; mechanical breaking does not break this bond. Therefore, within a crop rotation cycle, the earlier deep tillage is implemented, the more straw-derived carbon remains in the aggregate-encapsulated state, and the greater the amount released by disturbance. Conversely, the later the deep tillage is implemented, although the vulnerable portions have already fallen back, very little convertible carbon remains to be transported to lower layers and continue forming mineral-bound states. There is a suitable period in the middle of the crop rotation cycle, which is defined by two measurable quantities in this invention.

[0018] See Figure 1 This invention does not uniformly mix straw into the tillage layer. Instead, it uses a deep application and returning-to-field machine equipped with a double-layer fixed-depth furrow opener to create strips of two different soil depths within the same field. The strips located in the 10cm to 15cm soil layer are designated as shallow strips, and the strips located in the 25cm to 30cm soil layer are designated as deep strips. Only one strip is filled at each furrow opening location. Figure 1 The cross-section shows four consecutive strips arranged in alternating shallow, deep, shallow, and deep positions, with an adjacent strip spacing of 0.47 m, a strip width of 0.09 m, and a strip thickness of 0.025 m. Figure 1 The scale on the left indicates the soil depth. The dashed lines in the figure, which are 5cm apart, mark the boundaries of subsequent strata for sampling. Sampling extends from the surface down to 45cm. Figure 1 Two inverted triangle markers located in the 10cm to 15cm and 25cm to 30cm soil layers are designated as hydrothermal time measurement points, used to collect the temperature and moisture content of the soil layers where the two strips are located, respectively. Figure 1 The arc-shaped dashed arrow in the figure represents the vertical redistribution of straw carbon caused by deep tillage operations, which is the process that this invention aims to utilize at the appropriate time. Figure 1 The marked conditioning material is filled with strips filled with intersecting lines, while the unmarked conditioning material is filled with strips filled with dots.

[0019] Example 1 used alluvial soil plots in the North China Plain, with a light loam texture and a previous crop of maize. Before the experiment, undisturbed soil columns ranging from 0cm to 60cm were collected from each plot. Soil bulk density, cumulative soil mass, percentage of particles smaller than 0.02mm, organic carbon content in each carbon pool, and percentage of water-stable aggregates larger than 0.25mm were measured in 5cm layers. These values ​​formed the baseline. In this example, the bulk density of the 0cm to 20cm soil layer was 1.34g / cm³, and the bulk density of the 20cm to 45cm soil layer was 1.46g / cm³. 3 The mass percentage of particles smaller than 0.02 mm was 28.5%, and the mass percentage of water-stable aggregates larger than 0.25 mm was 38.6%. The sampling depth was 60 cm, while the calculation depth was 45 cm. A downward buffer was left when redefining the layer boundaries according to the equal mass soil layer method to avoid the need to push the layer boundaries to the unsampled soil layer due to changes in bulk density caused by deep tillage operations.

[0020] The conditioning material is prepared from corn stalks. After harvesting, the stalks are dried at 65℃ to constant weight and then graded into two particle size moduli using a hammer mill equipped with a two-stage sieve. The portion with a particle size of 0.5mm to 2mm is designated as the fine particle size moduli, and the portion with a length of 20mm to 50mm is designated as the coarse particle size moduli. A 2mm to 20mm gap is maintained between the two moduli; material falling into this gap is returned to the mill for reprocessing. This gap ensures that the material's state at the time of loading strictly corresponds to the definition of subsequent carbon stock operation: the fine particle size moduli can pass through a 2mm sieve and has a density of 1.8g / cm³. 3 The particles float in the heavy liquid and thus fall completely into the free particulate organic carbon pool at the time of filling; the coarse particle size modulus is retained by the 2mm sieve and thus falls completely into the coarse residual carbon pool at the time of filling. If a continuous particle size distribution is adopted, materials close to the 2mm boundary will be classified into different carbon pools during sieving due to differences in moisture content and operating force. The initial values ​​of the kinetic model can only be given by estimation, and the resulting deviation will be directly transmitted to each rate constant.

[0021] Labeled and unlabeled processed maize were obtained by parallel cultivation of the same maize variety in two separate cultivation chambers. The light, temperature, moisture, nutrient, and carbon dioxide concentration settings were identical in both chambers, differing only in the carbon-13 atomic fraction of the supplied carbon dioxide. At least three cultivation batches were produced, and the carbon-13 atomic fraction of the supplied carbon dioxide was exchanged between adjacent batches to offset any residual minor environmental differences between the chambers and avoid mistaking the chamber effect for an isotope effect. After harvest, the organic carbon mass fraction, carbon-nitrogen ratio, lignin mass fraction, cellulose mass fraction, hemicellulose mass fraction, and water-soluble organic carbon mass fraction of the labeled and unlabeled processed maize were measured in both coarse and fine particle size modalities. The relative deviation of each item under the same particle size modal was required to be no greater than 5%. The measured results of this embodiment are: organic carbon mass fraction 0.42, carbon-nitrogen ratio 63, lignin mass fraction 0.081, and water-soluble organic carbon mass fraction 0.052. The relative deviations of each item are all within 3%, which indicates that the two materials are only different in terms of carbon 13 atomic fraction.

[0022] Isotope specifications are based on the carbon-13 atomic fraction, and the measured values ​​are provided by an elemental analyzer coupled with a stable isotope ratio mass spectrometer. The instrument directly reports the deviation value relative to the international standard in parts per thousand (ppm), which needs to be converted to an isotope ratio value first. Then, the isotope ratios are converted into the carbon-13 atomic fraction. In the formula, The abundance ratio of carbon-13 to carbon-12 in the sample is given by a dimensionless quantity. The isotype ratio of the international reference material is taken as the valid international recommended value at the time of application; This is the deviation value reported by the instrument, expressed in parts per thousand. The carbon-13 atomic fraction is a dimensionless quantity. In this embodiment, the carbon-13 atomic fraction of the labeled prepared material is 0.0172, and the carbon-13 atomic fraction of the unlabeled prepared material is 0.01085.

[0023] The source tracing calculation is based on the carbon-13 atomic fraction instead of directly using the per mille deviation value because the relationship between the per mille deviation value and the mixing ratio is not linear. Within the natural abundance range, the deviation between the two is very small and can be approximated as linear. However, the carbon-13 atomic fraction of the labeled conditioning material is about 1.6 times that of the natural abundance. If linear mixing calculations are still performed using the per mille deviation value, the resulting labeled source carbon ratio will have a systematic deviation of several percentage points. After switching to the carbon-13 atomic fraction, the mixing relationship is strictly linear. Furthermore, since the labeled source carbon ratio is the ratio of the difference between two atomic fractions to the difference between two other atomic fractions, small variations in the international standard ratio are transmitted in the same direction in the numerator and denominator and cancel each other out. The impact on the calculation result is on the order of one ten-thousandth and can be ignored.

[0024] The strip application parameters are determined by back-calculation of the amount of carbon added per unit area. First, the carbon loading per unit strip length is calculated based on the strip cross-section and loading status. Then, the spacing between adjacent strips is calculated based on the amount of carbon input per unit area. In the formula, Carbon loading per unit strip length, in kg / m; This refers to the strip width, in meters (m). The thickness of the strip is expressed in meters (m). This refers to the bulk density of the filling material, expressed in kg / m³. 3 ; The organic carbon mass fraction of the conditioning material is a dimensionless quantity. The spacing between adjacent stripes, in meters; Carbon input per unit area, expressed in kg / hm² 2 In this embodiment, the strip width is 0.09m, the strip thickness is 0.025m, and the filling density is 150kg / m³. 3 With an organic carbon mass fraction of 0.42, the calculated carbon loading per unit strip length is 0.142 kg / m; the carbon input per unit area is 3000 kg / hm². 2 Reverse calculation yields an adjacent strip spacing of 0.47 m. Since the organic carbon mass fraction varies slightly depending on the coarseness-to-fineness ratio of the conditioning material, the adjacent strip spacing, strip width, and strip thickness are fixed during implementation, only at 140 kg / m². 3 Up to 160kg / m 3 Within the interval, the packing density is finely adjusted so that the carbon packing amount per unit strip length in each treatment is between 0.135 kg / m and 0.150 kg / m, so that the gradation effect will not be mixed with the change in strip spacing.

[0025] The depth of carbon return to the field is determined by the ratio of shallow to deep strips in every four consecutive strips, with three levels: 3:1, 2:2, or 1:3. Since only one strip is filled at each trench location, and the cross-section and bulk density of each strip are identical, the total number of strips, the amount of carbon input, and the trenching disturbance are completely consistent across the three configurations. The only difference lies in the depth distribution of the carbon input. If two strips are stacked vertically at the same location and their respective filling amounts are adjusted to achieve the depth distribution, the cross-section or bulk density of the strips in different configurations will inevitably change, and the effect of depth distribution will become intertwined with the effect of material compaction.

[0026] See Figure 2 Each cell is divided into shallow-marked subregions, deep-marked subregions, reference subregions, double-marked subregions, and blank subregions. Figure 2The comparison on the right shows the filling method for each subzone: In the shallow-marked subzone, only the shallow strips are filled with marked conditioner, while the deep strips are filled with unmarked conditioner; the deep-marked subzone is the opposite; in the reference subzone, both the shallow and deep strips are filled with unmarked conditioner; in the double-marked subzone, both strips are filled with marked conditioner; in the blank subzone, trenches are dug according to the same strip ratio and backfilled with an equal volume of native soil. The conditioner mass ratio, strip ratio, strip width and thickness, filling bulk density, carbon input, and all field operations are identical across the five subzones. The only difference between subzones is the strip position of the C13 label. Subzones are randomly arranged within the plot, with a buffer zone at least one machine working width between adjacent subzones to prevent tillage machinery from dragging marked conditioner from one subzone into an adjacent subzone.

[0027] Figure 2 The left side shows a magnified view of a shallowly marked sub-region. The marked conditioning material is not laid along the entire strip, but rather in 0.20m long marked sections at intervals, with a minimum distance of 0.50m between adjacent marked sections. Each marked section has a 0.15m long unmarked conditioning material protection section at both ends, and the remaining area is filled with unmarked conditioning material. During sampling, a section of soil 0.50m long, 0.50m wide, and 0.45m deep is excavated at the marked section. There is a definite relationship between the marked section length, the protection section length, and the sampling excavation length: before implementation, a preliminary test is conducted on the same field using an inert tracer to measure the maximum lateral migration distance under the adopted tillage method, which in this example is 0.15m; the protection section length is taken as this maximum lateral migration distance, and the sampling excavation length is the sum of the marked section length and twice the maximum lateral migration distance, which is exactly 0.50m. With this arrangement, carbon that migrates from the marked section will remain within the excavation area regardless of which side it moves to, and will not be misjudged as mineralization or leaching loss; at the same time, the marked sections on adjacent strips are staggered along the strip direction, and the marked conditioning material from another strip will not be mixed into one excavation window.

[0028] The amount of labeling and conditioning material used was calculated according to the above-described layout. The material mass of each label segment was the product of the strip cross-section, the label segment length, and the bulk density, approximately 0.068 kg. Nine label segments were set up for each type of strip in each subregion, corresponding to nine sampling time points. Phase 1 consisted of 27 plots: 27 shallow-position labeling subregions and 27 deep-position labeling subregions. Nine dual-position labeling subregions were set up in the first replicate of each treatment. Based on this, the required amount of labeling and conditioning material for Phase 1 was approximately 43.7 kg. In Phase 2, one additional set of paired label segments was set up in each plot for pre-plowing measurements, requiring approximately 5.5 kg, for a total of approximately 49 kg for the entire experiment. The labeling and conditioning material was only used for the label segments and not spread across the entire strip, keeping the tracer cost within an acceptable range. Subsequent kinetic and phase measurements were normalized using the label carbon input of the label segment itself; localized layout did not affect the comparability of the results.

[0029] Phase 1 was conducted under no-till conditions, with samples taken on days 0, 7, 30, 90, 180, 360, 720, 1080, and 1440 after conditioning. No basal tillage was performed during the tracking period, thus making the identified conversion rate constant a quantity determined by the conditioning, soil, and climate, independent of any planned deep tillage. If deep tillage were performed during the tracking period, the trajectory used to determine the timing of deep tillage would already include its influence, leading to self-referential criteria. The excavated soil was divided into nine 5cm layers, with the layer boundaries perfectly consistent between the entire treatment and the blank sub-region. Strips were assigned to their corresponding layers based on their absolute depth.

[0030] See Figure 3 Soil samples from each layer were separated into five non-overlapping carbon libraries as follows: Soil samples were first sieved through a 2mm sieve. Identifiable material residues on the sieve were manually removed, rinsed with deionized water, dried, weighed, and directly measured to obtain the coarse residue carbon library. The soil sample below the sieve was added to a sodium polytungstate solution with a density of 1.8 g / cm³ and centrifuged at 1000 times gravity for 15 min. The floating matter was the free particulate organic carbon library. The sediment was washed to remove the heavy liquid and then ultrasonically dispersed at an energy density of 300 J / mL. This density separation was repeated once, and the second floating matter was the closed-state particulate organic carbon library. The remaining heavy components were separated into fine-particle heavy-particle organic carbon libraries and coarse-particle heavy-particle organic carbon libraries at a 0.02mm boundary. The mass recovery rate of the five carbon libraries was required to be greater than 95%, and the carbon recovery rate to be between 95% and 105%. Samples that did not meet the standards were re-separated.

[0031] Figure 3Using the particle size of organic carbon and its resistance to mineralization as two conceptual coordinates, the transformation pathways among the five carbon pools mentioned above are illustrated. The coarse residual carbon pool is crushed into the free particulate organic carbon pool; part of the free particulate organic carbon pool is coated by aggregates into the closed-state particulate organic carbon pool, part is directly absorbed into the fine-grained heavy component organic carbon pool through microbial assimilation and mineral surface adsorption, and another part combines with the surface of sand grains to enter the coarse-grained heavy component organic carbon pool; the closed-state particulate organic carbon pool and the coarse-grained heavy component organic carbon pool can further transform into the fine-grained particulate organic carbon pool; the coarse residual carbon pool, the free particulate organic carbon pool, and the closed-state particulate organic carbon pool are also lost through mineralization and leaching, respectively. Figure 3 The diagonal thick dashed line in the figure divides the five carbon pools into tillage-sensitive zones and tillage-insensitive zones. The closed-state particulate organic carbon pool is located on one side of the sensitive zone, and the arrow of deep tillage disturbance points directly to this carbon pool; the fine-grained heavy component organic carbon pool and the coarse-grained heavy component organic carbon pool are located on one side of the insensitive zone. Figure 3 The text also indicates which carbon library each of the two quantities subsequently used as phase coordinates corresponds to.

[0032] The heavy components obtained from density separation are not all organic carbon bound to the surface of fine-grained minerals; some are adsorbed onto the surface of sand grains. Therefore, the heavy components are further divided into two carbon libraries based on a thickness of 0.02 mm. The organic carbon library of the fine-grained heavy components corresponds to the stabilization mechanism of mineral surface binding. Its stability is independent of the integrity of the aggregates and is therefore insensitive to mechanical disturbances. The organic carbon library of the coarse-grained heavy components is analyzed separately, which ensures the mass closure of the five carbon libraries and avoids overestimating the degree of stabilization by including sand-bound carbon in the mineral-bound proportion. The carbon library of the coarse residues is directly measured without removing inorganic carbon because plant residues themselves do not contain carbonates, and acid washing would instead leach out water-soluble organic carbon and change the composition. For soil samples containing carbonates, inorganic carbon is removed by hydrochloric acid vapor fumigation, as soaking and acid washing would also result in the loss of soluble organic carbon.

[0033] After processing, the organic carbon content and per mille deviation of each carbon library were measured and converted into carbon-13 atomic fraction. Then, the carbon proportion of the labeled source was calculated according to the mixing relationship between the two endmembers. In the formula, To indicate the proportion of carbon from the source, it is a dimensionless quantity; The carbon-13 atom fraction of the carbon library being tested; The carbon 13 atom fraction in the reference subregion for the same cell, the same sampling time point, the same layer, and the same carbon library; This refers to the carbon 13 atomic fraction of the labeled conditioning material. When calculating the source of shallow bands, the carbon library is taken from the shallow labeled subregion; when calculating the source of deep bands, it is taken from the deep labeled subregion. The labeled source carbon reserves are then calculated using the following formula. In the formula, To indicate the source carbon reserves, the unit is g / m³. 2 ; The value represents the mass content of organic carbon in the whole soil layer of this carbon pool, expressed in g / kg. This is the unit weight of the soil layer at this depth, expressed in g / cm³. 3 ; The layer thickness is taken as 5 cm; the coefficient 10 is derived from unit conversion. This formula is not applicable to coarse residue carbon pools; their reserves are calculated by dividing the product of the dried mass of the recovered residue and the mass fraction of carbon in the residue by the area of ​​the sampling unit. The marked source carbon reserves of each carbon pool are summed over the entire profile from 0 cm to 45 cm and normalized by the marked carbon input amount of that marked segment. After normalization, they are recorded sequentially as follows: to .

[0034] Background subtraction uses the carbon pool of the same name in the reference subregion, rather than the blank subregion. The reference subregion is completely identical to the shallow-marked subregion and the deep-marked subregion in terms of material chemical composition, particle size distribution, band number ratio, carbon input, band geometry, backfilling, and all field operations. The only difference is the carbon-13 atomic fraction of the loaded material. After subtracting the two, four items are simultaneously subtracted: background soil organic carbon, root carbon of the current crop and rhizosphere deposited carbon, vertical redistribution of background carbon caused by tillage operations, and the primary organic carbon activation effect caused by the input of the same amount of organic material itself. If a blank subregion without any material is used as the background, the first three items can be subtracted, but the activation effect cannot be subtracted, and the resulting marked carbon stock will be systematically overestimated.

[0035] The consistency of the results was confirmed by two checks. The first check was the additivity check, which required that the relative deviation between the labeled carbon reserves obtained from the two-position labeled subregion and the sum of the reserves obtained from the shallow-position labeled subregion and the deep-position labeled subregion should not exceed 10%. In this embodiment, the measured value was 6.3%, indicating that there was no significant mutual interference between the destinations of the input carbon from the two types of bands, and they could be obtained separately and then synthesized. The second check was the diagnostic check, which utilized the natural abundance difference between corn stalks, which are C4 crop residues, and the tested soil, which has long been dominated by C3 crops. The carbon 13 atomic fraction was calculated independently by dividing the difference between the carbon 13 atomic fraction of the same carbon library in the reference subregion and the blank subregion by the difference between the carbon 13 atomic fraction of the same carbon library in the unlabeled conditioning material and the blank subregion. This channel is systematically higher due to the influence of root carbon and excitation effects. In this embodiment, it is 11.8% higher than the labeled channel. It is only used for diagnosis of order of magnitude consistency and is not used as a criterion for data acceptance.

[0036] The kinetic identification uses hydrothermal time, rather than calendar time, as the independent variable. Hydrothermal time is accumulated according to the following formula. In the formula, Hydrothermal time, expressed in °C·d; The average daily temperature of the soil layer on that day, in °C; As an indicative variable, the soil moisture content of the soil layer on a given day is taken as 1 if it is not less than 50% of the field capacity, and as 0 otherwise; the summation iterates through each day from the date of filling the conditioning material. The summation is based on the traversed day sequence. The shallow strip uses temperature and moisture content data from the 10cm to 15cm soil layer, and the deep strip uses data from the 25cm to 30cm soil layer. Figure 1 Data were collected from two hydrothermal time measurement points. The annual temperature amplitude and water content of the soil layers where the two types of strips are located are not the same, and therefore the biologically effective time corresponding to the same calendar day is also different. Therefore, hydrothermal time is accumulated separately according to the soil layer depth of each strip. By using hydrothermal time, the same set of rate constants can simultaneously describe the transformation process of the two types of strips, and the difference between the two types of strips is attributed to the difference in the hydrothermal time accumulation rate. In this embodiment, the annual accumulation of the shallow strip is approximately 2400℃·d, and that of the deep strip is approximately 2000℃·d.

[0037] See Figure 4 For each type of stripe, the rate constant is identified according to the following first-order kinetic equations: ; ; ; ; In the formula, to The labeled carbon reserves in the coarse residual carbon library, the free particulate organic carbon library, the closed storage particulate organic carbon library, the fine heavy component organic carbon library, and the coarse heavy component organic carbon library have all been normalized to dimensionless quantities based on the labeled carbon input amount. The conversion rate constant from the coarse residual carbon library to the free particulate organic carbon library; , and The conversion rate constants, in order, represent the transformation from the free particulate organic carbon library to the closed-state particulate organic carbon library, the fine-particle heavy component organic carbon library, and the coarse-particle heavy component organic carbon library. is the conversion rate constant from the closed-state particulate organic carbon library to the fine-particle heavy component organic carbon library; The conversion rate constant from the coarse-grained heavy component organic carbon library to the fine-grained heavy component organic carbon library; , , , and The following are the rate constants for the loss of the corresponding carbon pools in the form of mineralization and leaching, respectively; the units of the above rate constants are all reciprocals of °C·d.

[0038] The initial values ​​of the equations are directly given by the material composition: the initial value of the coarse residual carbon pool is taken as the carbon content of the coarse particle size mode filled in the marked segment; the initial value of the free particulate organic carbon pool is taken as the carbon content of the fine particle size mode filled in the marked segment; and the initial values ​​of the other three carbon pools are 0. Sampling on day 0 is used to verify whether this correspondence holds. In this embodiment, the measured reserves of the coarse residual carbon pool and the free particulate organic carbon pool on day 0 are 0.98 times and 0.97 times the theoretical filling amount, respectively. The other three carbon pools are below the detection limit, indicating that the bimodal particle size composition indeed allows the initial values ​​to be given by weighing without estimation. The fitting uses a nonlinear mixed-effect model with the natural logarithm of each rate constant as a parameter. The ratio of coarse to fine mass and the ratio of the number of strips in the conditioning material act on the fixed effects on a logarithmic scale, while the small plots are random effects. Parameterizing on a logarithmic scale allows the rate constants to naturally take positive values ​​without additional boundary constraints and makes the effect of the treatment factors on the rate appear as a multiple relationship, facilitating comparisons across soil types.

[0039] Before fitting, a parameter identifiability test is performed, which involves calculating the rank and condition number of the sensitivity matrix. Rank-deficient parameters are fixed one by one according to the Akaike information content criterion before refitting, and the sensitivity analysis results are reported. In this embodiment... and The sensitivity within the observation window was significantly low. After comparing two models—one containing this parameter and the other with it fixed at 0—according to the Akaike information content criterion, the model was retained. And The value is fixed at 0. The final fitted main rate constant is: 0.00015 0.00015 0.00008 0.00009 The value is 0.00030, and the units are all the reciprocals of ℃·d. The relative width of the 95% confidence interval obtained by self-sampling 1000 times with the complete cell as the unit is within 25%.

[0040] Figure 4 The graph shows the variation of labeled carbon reserves of five carbon pools over calendar time at a coarse-to-fine mass ratio of 1:1, where (a) represents shallow bands and (b) represents deep bands. The horizontal axis of both graphs represents calendar time in days; the vertical axis represents normalized labeled carbon reserves, which is dimensionless. Figure 4 In the study, the coarse residual carbon pool and the free particulate organic carbon pool continuously decreased from the date of loading, while the fine heavy component organic carbon pool increased monotonically. The closed-system particulate organic carbon pool first increased and then decreased, reaching its peak around day 550 in (a) and delayed to around day 660 in (b) due to slower hydrothermal accumulation. This pattern of initial increase followed by decrease in the closed-system particulate organic carbon pool makes it suitable for describing the costs of deep tillage operations.

[0041] The normalized trajectories of the two stripes were weighted and synthesized using calendar time as a common index. In the formula, This is calendar time, in days (d). For the synthesized first Normalized labeled source carbon storage of each carbon pool; and They are respectively the superficial bands and the deep bands. Normalized labeled source carbon storage for each carbon pool, with each having its own hydrothermal time as the independent variable; and These are the cumulative functions of hydrothermal time with calendar time for shallow and deep stripes, respectively; and The ratios of carbon input carried by the shallow and deep strips are respectively, and their sum is 1. The hydrothermal time of the two strips advances at different rates, and the same hydrothermal time value corresponds to different calendar dates on the two strips. If the hydrothermal time is directly added as an index, the resulting curve does not correspond to any real time. Therefore, the synthesis must first be converted to calendar time.

[0042] Based on the synthetic trajectory, three quantities are defined: mineral binding recovery rate. Closed-loop protection recovery rate Overall retention rate In the formula, Mineral binding recovery rate, which is the ratio of labeled source carbon reserves to labeled carbon input in the organic carbon library of fine-grained heavy components; The closed-loop recovery rate is the ratio of labeled carbon reserves to labeled carbon input in the closed-loop particulate organic carbon pool. The total retention rate is the ratio of the sum of the carbon reserves from the five carbon pools to the amount of labeled carbon input; all three are dimensionless. The mineral-bound recovery rate represents the benefit side of deep tillage operations; the higher the value, the more straw-source carbon has entered a state of insensitivity to tillage. The closed-cell recovery rate represents the cost side of deep tillage operations; the higher the value, the more straw-source carbon is currently protected by the integrity of the aggregates, and the greater the amount released once deep tillage is implemented. The total retention rate is used as a constraint to prevent deep tillage operations from being postponed until a large amount of straw-source carbon has been lost.

[0043] See Figure 5 The horizontal axis represents the mineral binding recovery rate, and the vertical axis represents the closed-loop protection recovery rate. Both are dimensionless. The three curves in the figure correspond to the phase coordinate trajectories of conditioning materials with coarse-to-fine mass ratios of 3:1, 1:1, and 1:3, respectively, when the strip number ratio is 2:2. The arrows on the curves indicate the direction of time progression. Because the mineral binding recovery rate increases monotonically while the closed-loop protection recovery rate first increases and then decreases, the trajectory in this plane shows a shape that first moves to the upper right, crosses the inflection point, and then turns to the lower right. Figure 5The circular, triangular, and square markers indicate the locations of the basic tillage operation windows for the 1st, 2nd, and 3rd years, respectively. The star-shaped marker represents the calibrated phase coordinate target value. The two dashed lines and the shaded area together define the feasible region, which is the area where the mineral binding recovery rate is not lower than the target value and the closed-loop protection recovery rate is not higher than the target value. Deep tillage operations should be arranged in the first basic tillage operation after the trajectory first enters the feasible region. Figure 5 It is evident that only the third-year marker of the trajectory with a thickness-to-weight ratio of 1:1 falls within the feasible region.

[0044] The combination of mineral-bound recovery rate and closure-protected recovery rate was selected instead of the relative proportions of fine-grained heavy component organic carbon pools in the soil matrix, after coordinate identifiability calculations. If relative proportions are paired with mineral-bound recovery rate, they share the same molecule, resulting in a correlation coefficient of 0.83 across all design points, a variance inflation factor of 3.19, a design matrix condition number of 31.0, a relatively small convex hull area spanned by the measured coordinates, and a minimum detectable effect of 294 kg / hm² for the quadratic and interaction terms of the response surface. 2 In terms of carbon dioxide equivalent, this exceeds the potential inter-treatment variability, and the response surface effectively degenerates into a one-dimensional curve. After pairing the mineral-bound recovery rate with the closed-loop protection recovery rate, the correlation coefficient decreased to 0.42, the variance inflation factor was 1.21, the convex hull area increased by approximately 60%, and the minimum detectable effect of the interaction term decreased to 82 kg / hm². 2 The two-dimensional response surface can be truly identified. This is because the recovery rate of the closed-loop protection first increases and then decreases over time, and there is no monotonic correspondence between it and the monotonically increasing mineral-bound recovery rate. Therefore, it can provide a second independent direction for the response surface.

[0045] Phase 2 involved a fully crossover experiment with three factors: the coarse-to-fine ratio of conditioning materials, the proportion of strips, and the tillage sequence. A total of 27 treatments were implemented, replicated three times, and initiated in three consecutive years, comprising 81 plots. In each of the four basal tillage operations within a single crop rotation cycle, the three tillage sequences included only one deep tillage operation, scheduled for years 1, 2, and 3 respectively. The remaining operations were rotary tillage, ensuring that the only difference between the three levels was the year in which the deep tillage operation occurred. This year-by-year initiation provided treatment-independent interannual climate variability, thus the measured phase coordinates were no longer deterministic functions of the treatment labels, and the response surface could be identified simultaneously with the main effects of the treatments. Furthermore, year-by-year initiation provided conditions for external validation, allowing the results of the first-initiated block to predict the results of the subsequent-initiated blocks, thereby predicting the robustness of the residual test calibration. When the field area is limited, the D-optimal design, which can estimate all two-factor interaction terms and has no less than 24 different design points, can be used instead. The number of design points must be no less than the number of parameters to be estimated in the model and allowance for fit testing.

[0046] In Phase Two, each plot had an additional pairing marker section during the filling of conditioning materials. One marker section was filled with labeled conditioning materials, while the other, paired with unlabeled conditioning materials, served as a reference. The materials in both marker sections were in direct contact with the soil. Excavation was conducted within 7 days prior to deep tillage in each plot, and the materials were graded as described above, with residues and surrounding soil collected simultaneously. This allowed for direct measurement of the mineral binding recovery rate and the carbon sequestration recovery rate at the time of deep tillage. Phase Two involved re-measurement instead of reading from the Phase One trajectory because the number of rotary tillage operations experienced by the three tillage sequences before deep tillage varied. The plot deep-tilled in the first year had not undergone prior rotary tillage, while the plots deep-tilled in the second and third years had undergone one and two rotary tillage operations, respectively, resulting in a carbon pool state that deviated from the trajectory under no-till conditions. Using direct measurement, the phase coordinates aligned with the actual history of each plot. If mesh bags are used for isolation, the contact between materials and minerals, the entry of soil animals, and the transfer of carbon to the soil outside the bag will be altered, and the resulting trajectory cannot represent the actual process within the strip.

[0047] Net carbon sequestration was calculated using the difference-in-differences and equal-mass soil layer method. In the formula, Net carbon sequestration, expressed in kg / hm², in carbon dioxide equivalent. and These represent the organic carbon storage in the 0cm to 45cm depth of the treated plot at the accounting time point and the baseline time point, respectively. and These represent the same-named reserves of blank sub-regions within the same community at the corresponding time points, all in kg / hm², expressed as carbon; the coefficient 44 to 12 is the mass ratio of carbon to carbon dioxide. The carbon dioxide equivalent of the cumulative increase in methane emissions, converted at 27%. The carbon dioxide equivalent of the cumulative increase in nitrous oxide emissions, converted to 273. The emissions are calculated based on an increase in diesel consumption for tillage operations, converted at 2.68 kg of carbon dioxide per liter; all items are calculated based on the same area and the same time period. The baseline soil mass is the minimum cumulative soil mass from 0cm to 45cm in all plots. Taking the minimum value ensures that all plots only need to be contracted upwards when redefining the strata boundaries, without having to be pushed downwards.

[0048] Soil carbon dioxide flux is not deducted in the above formula. The change in organic carbon storage is already a comprehensive result of carbon input, decomposition, and release during that period. If the sequestration is calculated based on the increase in storage first and then subtracted from the soil carbon dioxide flux, the same decomposition loss will be calculated twice. Soil carbon dioxide flux is separately statistically analyzed in this invention to verify whether the material mineralization rate and activation process are consistent with the kinetic identification results, and it is not included in the net carbon sequestration. Methane and nitrous oxide are different; they do not constitute part of the soil organic carbon storage, and their emission increments must be deducted separately. Greenhouse gas fluxes are measured every 7 to 10 days throughout the growing season using a static chamber coupled with gas chromatography, and are measured more frequently for 5 consecutive days after tillage, fertilization, daily precipitation greater than 10 mm, irrigation, and freeze-thaw events to avoid missing emission peaks caused by fertilization and precipitation due to short-term monitoring after tillage.

[0049] See Figure 6 Using net carbon fixed as the dependent variable, and the primary and secondary terms of mineral-bound recovery rate and closed-loop protection recovery rate, as well as their interaction term, the main effects of the coarse-to-fine mass ratio of conditioning materials and the ratio of strip number, as well as their interaction term, carbon loading per unit strip length, and the random intercept of the start-up year as independent variables, a response surface model was fitted. The cultivation sequence was not included as a factor in the response surface model. This invention argues that the effect of the cultivation sequence is entirely transmitted through phase coordinates. To verify this premise, an alternative model was fitted with the cultivation sequence as a factor but without phase coordinates, and the two were compared according to the Akaike information criterion and cross-validation. In this embodiment, the Akaike information criterion value of the response surface model containing phase coordinates is 18.4 lower than that of the alternative model, and the root mean square error of leave-one cross-validation is 23% lower, indicating that the effect of the cultivation sequence can indeed be summarized by phase coordinates.

[0050] Before fitting, a coordinate identifiability test is performed, requiring that the variance inflation factors of both the mineral binding recovery rate and the closed-loop protection recovery rate are less than 2, the design matrix be of full rank and the condition number less than 30, and the constraint optimum fall within the convex hull formed by the measured phase coordinates. In this embodiment, the variance inflation factor is 1.21, the design matrix has a rank of 14 and is of full rank, the condition number is 25.1, and the constraint optimum is located inside the convex hull; all indicators meet the requirements. Figure 6 The horizontal axis represents the mineral binding recovery rate, and the vertical axis represents the closed-loop protection recovery rate; both are dimensionless. The closed curve is the net carbon solidification contour line, with the unit of the labeled values ​​being kg / hm², expressed as carbon dioxide equivalent. The scatter plots represent the measured phase coordinates of the 81 experimental plots. The thick dashed line represents the convex hull of the measured phase coordinates. The star-shaped markers and their error bars represent the constrained optimal points and their 95% confidence intervals.

[0051] The constraints are: the lower confidence limit of the total retention rate is not lower than a pre-set threshold, and the mass ratio of water-stable aggregates is higher than the baseline value. In this embodiment, the threshold is determined to be 0.35 by the baseline and pre-experiment, which is used to limit deep tillage operations to a period when there is still sufficient convertible carbon from straw. Under this constraint, the constrained optimal points of the response surface model are obtained as follows: mineral-bound recovery rate of 0.130 and closed-loop protection recovery rate of 0.078, with 95% confidence intervals of 0.120 to 0.140 and 0.071 to 0.085, respectively. The target values ​​of the phase coordinates are taken from the conservative side of this confidence interval, that is, the target mineral-bound recovery rate is taken from the lower confidence limit of 0.120, and the target closed-loop protection recovery rate is taken from the upper confidence limit of 0.085. Taking the conservative side can avoid the following situation: field measurements inevitably contain errors. If the target value is taken at the optimal point itself, when the measured value fluctuates around the optimal point, the deep tillage operation time will jump back and forth between two adjacent basic tillage operations.

[0052] Based on the aforementioned target values, the mineral binding recovery rate reached 0.120 on day 816 after loading, and the closed-loop protection recovery rate dropped to 0.085 on day 940. The earliest time when both targets were met simultaneously was day 940. The four basic tillage operations occurred on days 365, 730, 1095, and 1460, respectively. Therefore, the first feasible basic tillage operation window was the third year, i.e., the tillage sequence was rotary tillage, rotary tillage, deep tillage, and rotary tillage. In actual implementation, the above calculations were not relied upon. Instead, paired marker sections were excavated and measured directly before each candidate basic tillage operation window. The window that first simultaneously met both target values ​​was designated as the deep tillage operation window, and all other basic tillage operations used rotary tillage. If the two target values ​​cannot be met simultaneously within a single crop rotation cycle, then all basic tillage operations in that cycle will be carried out using rotary tillage, and deep tillage operations will be postponed to the next crop rotation cycle, rather than being forced to carry out deep tillage operations in the final window.

[0053] The number of deep-seated bands is selected based on the highest net carbon fixation predicted by the response surface methodology, under the constraint that the organic carbon saturation deficit of fine-grained heavy components is not less than 0. The saturation deficit is calculated using the following formula. In the formula, For saturation deficiency; This is an estimated capacity value for the organic carbon in the fine-grained heavy components of this layer; This refers to the measured organic carbon content of fine-grained heavy components in this soil layer; all three values ​​are in g / kg, expressed as the mass of carbon contained in one kilogram of soil. The capacity estimate is given by quantile regression from a regional broad-grained gradient reference database, or it can be determined using direct mineral adsorption experiments. In this example, the capacity estimate for the 25cm to 30cm soil layer is 14.6 g / kg, the measured organic carbon content of fine-grained heavy components is 9.8 g / kg, the saturation deficit is 4.8 g / kg, and its 95% confidence limit is 1.9 g / kg, thus the constraint is valid. This item is used as a propensity score constraint in the optimization process, rather than as a separate selection criterion, because the capacity estimate is derived from empirical relationships, and the reliability of its extrapolation for some soil types needs further verification.

[0054] The effects were calculated after implementing the above scheme for a complete 4-year crop rotation cycle. Under identical conditions, the net carbon sequestration for the three tillage sequences was as follows: 850 kg / hm² for deep tillage in the first year. 2 The yield was 1040 kg / hm² in the second year. 2 The rate was 1260 kg / hm² in the third year. 2 All figures are expressed in carbon dioxide equivalents. The second year's basal tillage operation window falls near the peak of the closed-loop conservation recovery rate. At this time, the amount of straw source carbon encapsulated in aggregates is the highest, and the amount released by deep tillage is also the largest. The increase in benefits is offset by the increase in costs, and its value is therefore not significantly different from the first year. If deep tillage is further postponed to the fourth year, the net carbon sequestration predicted by the response surface model falls back to 850 kg / hm², comparable to the level of the first year. At this point, the total retention rate has dropped to 0.325, below the set threshold of 0.35, and the constraint is not valid. This non-monotonic change in net carbon sequestration with the year of deep tillage indicates that the suitable period is within the rotation cycle rather than at its end; control measures based solely on the number of deep tillage cycles or the depth of deep tillage cannot achieve this effect.

[0055] The effect of particle size distribution was also demonstrated. Under the condition that deep tillage was consistently performed in the third year, the net carbon sequestration was 1110 kg / hm² when the coarse-to-fine particle size ratio was 3:1. 2 At a ratio of 1:1, the yield is 1260 kg / hm. 2 At a ratio of 1:3, the concentration is 1000 kg / hm². 2 All values ​​are expressed in carbon dioxide equivalent. When the proportion of coarse-grained mode is too high, the release of carbon from straw source is slow, and the mineral binding recovery rate is only 0.107 in the third year window, failing to reach the target value. When the proportion of fine-grained mode is too high, the release is too fast, and the closed-loop protection recovery rate is still 0.089 in the third year window, exceeding the target value. The interaction term between the coarse-fine mass ratio and the phase coordinate in the response surface model is significant, indicating that the two must be optimized together; simply superimposing the optimal values ​​of each does not yield the best result. When the particle size distribution is 1:3 and deep tillage is arranged in the second year, the net carbon sequestration drops to 95 kg / hm².2 The carbon fixation effect has essentially disappeared, further illustrating the cost of phase mismatch.

[0056] Regarding soil improvement, the proportion of water-stable aggregates with a particle size greater than 0.25 mm increased from 38.6% to 46.4% of the baseline. At the end of one crop rotation cycle, the labeled carbon storage was 900 kg / hm². 2 In terms of carbon, it accounts for 0.300% of the carbon input per unit area; the organic carbon storage calculated using the equal mass soil layer method from 0cm to 45cm is 848kg / hm² higher than that of the blank sub-region. 2 In carbon terms, the change in primary soil organic carbon, calculated from the difference between the net increase and the labeled source carbon storage, is -52 kg / hm². 2 This indicates the presence of a mild stimulating effect. The net increase, converted at a ratio of 44 to 12, is 3109 kg / hm². 2 The carbon dioxide equivalent, minus the cumulative increase in methane emissions of 6 kg / hm² 2 The cumulative increase in nitrous oxide emissions is 1820 kg / hm². 2 Increased diesel consumption for tillage operations by 23 kg / hm 2 After adjusting for carbon dioxide equivalent, the net carbon sequestration is 1260 kg / hm². 2 The increase in nitrous oxide emissions offsets most of the increase in organic carbon; if this is not taken into account, the carbon sequestration effect will be significantly overestimated. As a comparison, a conventional treatment was implemented on the same field: straw was pulverized and then completely buried in a 15cm rotary tillage layer, with the same carbon input of 3000 kg / hm². Deep tillage was consistently scheduled in the first year of the crop rotation cycle, resulting in a net carbon sequestration of 430 kg / hm². 2 The water-stable aggregates accounted for 41.2% of the total mass of carbon dioxide equivalent. Under the same conditions of total straw input and number of deep tillage operations, the net carbon sequestration of this invention is approximately 2.9 times that of the comparative example.

[0057] Example 2 illustrates the implementation of the invention on other soil types. The test plot was sandy loam black soil from the Huang-Huai region, with a medium loam texture. The percentage of particles smaller than 0.02 mm was 41.2%, and the estimated organic carbon capacity of the fine-grained heavy components in the 25-30 cm soil layer was 19.3 g / kg. This soil had a high content of fine-grained minerals and a significant saturation deficit. The optimization result was to increase the number of deep-seated bands to 3 out of every 4 consecutive bands, i.e., a band ratio of 1:3. The measured organic carbon mass fraction of the conditioning material was 0.404, from which the spacing between adjacent bands was calculated to be 0.45 m, while the geometric parameters of the other bands remained unchanged. The annual cumulative hydrothermal time of this field is slightly lower than that of Example 1, with the shallow stripe at approximately 2250℃·d and the deep stripe at approximately 1880℃·d. The calibrated phase coordinate target values ​​are 0.114 for mineral binding recovery and 0.091 for storage and conservation recovery, which are close to but not the same as those in Example 1. The earliest time when both target values ​​are simultaneously met is the 880th day, and the deep tillage operation window still falls in the 3rd year. The target values ​​need to be calibrated separately according to soil type. When used across soil types, the calibration should be repeated on no less than 3 fields with different textures or mineral compositions.

[0058] As an optional implementation method, when the implementing unit lacks the conditions for preparing carbon-13 labeled materials, only the natural abundance channel can be retained. This involves using unlabeled corn stalks as the sole material, with the corresponding carbon library in the blank subregion as the background, and calculating the source carbon ratio according to the same atomic fraction mixing relationship. In this case, the source tracing sensitivity is significantly reduced, and the excitation effect cannot be subtracted. The resulting phase coordinates have a systematically high bias. Therefore, it is necessary to first calibrate the bias using the labeled channel on the same soil type, and then correct the results of the natural abundance channel accordingly. This implementation method is suitable for already calibrated promotion fields, and its cost is far lower than a complete labeling scheme.

[0059] As an alternative implementation, hydrothermal timing data can be automatically collected by temperature and moisture sensors buried in soil layers of 10cm-15cm and 25cm-30cm. Data is accumulated daily by a field data logger, and an alert is issued when the accumulated value reaches a preset threshold, guiding sampling and tillage operations. Furthermore, the length of the coarse-grained mode can be adjusted within the range of 20mm-50mm, and the particle size of the fine-grained mode can be adjusted within the range of 0.5mm-2mm. As long as a 2mm-20mm particle size gap is maintained between the two modes, a one-to-one correspondence between the initial values ​​and the carbon pool operation definition can be maintained. The length of the crop rotation cycle can also be changed from 4 years to 3 or 6 years, in which case the number of basic tillage operations changes accordingly, while deep tillage operations are still arranged according to the window where both target values ​​are simultaneously met for the first time. For planting systems such as rice-wheat rotation with alternating flooding and drying, the moisture content criterion for hydrothermal timing needs to be changed to redox potential, while the other steps remain unchanged.

[0060] In this invention, the composition parameters of the conditioning material, the strip application parameters, and the tillage sequence constitute a unified output. The final output includes the bimodal particle size distribution and coarse-to-fine mass ratio of the conditioning material, the bulk density and carbon loading per unit strip length, the spacing between adjacent strips and the ratio of shallow to deep strips, the hydrothermal time threshold corresponding to the subsequent crop sowing sequence, the target phase coordinate value, and the tillage sequence arrangement within one crop rotation cycle determined by these parameters. In production, subsequent straw return and tillage operations are implemented according to this set of parameters, and the response surface model is updated with the current measured data at the end of each crop rotation cycle. When the 95% confidence lower limit of the net carbon fixation improvement of the predicted optimal combination relative to the currently implemented combination is greater than 0, the parameter combination is changed; otherwise, the current combination is maintained, and the iteration is performed for a maximum of three crop rotation cycles.

[0061] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention, such as replacing corn stalks with other C4 crop stalks, appropriately adjusting the soil depth of the two strips within the tillage and plow pan ranges, or replacing the nonlinear mixed-effect model with other parameter estimation methods with equivalent identification capabilities, should be included within the scope of protection of the present invention.

Claims

1. A method for soil improvement and carbon sequestration through synergistic regulation of straw return to the field and tillage sequence, characterized in that, Includes the following steps: Corn stalks were graded into fine and coarse particle size modes with non-overlapping particle size ranges and compounded into conditioning materials according to a set mass ratio. Labeled and unlabeled conditioning materials with different carbon 13 atom fractions but the same particle size mode composition and chemical composition were prepared. Within the same field, shallow and deep strips with different soil depths are cut according to the set number of strips and filled with conditioning material. At the same time, the field is divided into three sub-zones. The shallow marked sub-zone is filled with marked conditioning material only in the shallow strips, the deep marked sub-zone is filled with marked conditioning material only in the deep strips, and the reference sub-zone is filled with unmarked conditioning material. The application parameters for the three sub-zones are the same as those for field operations. Before deep tillage, soil samples were collected from each sub-region in layers to separate multiple non-overlapping carbon pools, including closed-state granular organic carbon pools and fine-grained heavy component organic carbon pools. Using the carbon 13 atomic fraction of the carbon pool of the same name in the reference sub-region as the baseline, the labeled source carbon reserves of each carbon pool were calculated based on the carbon 13 atomic fraction. The ratio of the labeled source carbon reserves of the fine-grained heavy component organic carbon pool and the closed-state granular organic carbon pool to the labeled carbon input was used as the mineral binding recovery rate and the closed-state protection recovery rate, respectively. The target value of the phase coordinate, which is composed of the mineral binding recovery rate and the closed storage protection recovery rate, was obtained by calibration using the response surface model. Deep tillage was scheduled to be carried out in the first basic tillage operation after the mineral binding recovery rate and the closed storage protection recovery rate reached the target value of the phase coordinate. Rotary tillage was used for the remaining basic tillage operations in the rotation cycle.

2. The method according to claim 1, characterized in that, The fine-particle size modulus has a particle size of 0.5 mm to 2 mm, and the coarse-particle size modulus has a segment length of 20 mm to 50 mm. A particle size gap of 2 mm to 20 mm is set between the two moduli. The mass ratio of the coarse-particle size modulus to the fine-particle size modulus is 3:1, 1:1, or 1:

3. The carbon-13 atomic fraction of the labeled conditioned material is 0.0165 to 0.0200, and the carbon-13 atomic fraction of the unlabeled conditioned material is 0.01080 to 0.01092. Under the same particle size modulus, the relative deviations between the labeled and unlabeled conditioned materials in terms of organic carbon mass fraction, carbon-nitrogen ratio, lignin mass fraction, cellulose mass fraction, hemicellulose mass fraction, and water-soluble organic carbon mass fraction are all no greater than 5%.

3. The method according to claim 1 or 2, characterized in that, Labeled and unlabeled processed materials were harvested from two cultivation chambers of the same maize variety, which were cultivated in parallel in two chambers with different carbon 13 atomic fractions of carbon dioxide but identical other cultivation conditions. There were no fewer than three cultivation batches, and the carbon 13 atomic fractions of carbon dioxide supplied by the two cultivation chambers were exchanged between adjacent batches.

4. The method according to claim 1, characterized in that, The shallow strips are located in the 10cm to 15cm soil layer, and the deep strips are located in the 25cm to 30cm soil layer. The spacing between adjacent strips is 0.47m, the strip width is 0.09m, the strip thickness is 0.025m, the bulk density is 140kg / m3 to 160kg / m3, and the carbon loading per unit strip length is 0.135kg / m to 0.150kg / m. The carbon input per unit area is 2900kg / hm2 to 3200kg / hm2. The ratio of the number of shallow strips to deep strips in every 4 consecutive strips is 3:1, 2:2, or 1:

3.

5. The method according to claim 1, characterized in that, The experiment was further divided into two sub-regions: a double-marked sub-region and a blank sub-region. In the double-marked sub-region, both shallow and deep strips were filled with marked conditioning material. In the blank sub-region, trenches were dug according to the same strip ratio and backfilled with an equal volume of original soil. The sub-regions were randomly arranged, with an isolation zone at least one machine operating width between adjacent sub-regions. The marked conditioning material was laid out along the strips in 0.20m long marked sections, directly contacting the soil. The distance between adjacent marked sections was at least 0.50m. Each marked section had a 0.15m long unmarked conditioning material protection section at both ends. The sampling excavation length was at least the sum of the marked section length and twice the maximum lateral migration distance, which was determined by a preliminary test using an inert tracer. Additivity verification was performed when the relative deviation between the carbon reserves from the two-position labeled subregion and the sum of the carbon reserves from the shallow-position labeled subregion and the deep-position labeled subregion was no greater than 10%. Diagnostic verification was performed when the carbon 13 atomic fraction of the reference subregion and the blank subregion corresponding carbon library was divided by the carbon 13 atomic fraction of the unlabeled conditioning material and the blank subregion corresponding carbon library.

6. The method according to claim 1, characterized in that, Carbon pool separation and source carbon storage calculation were performed as follows: Soil samples were sieved through a 2mm sieve. Identifiable residues on the sieve were picked out, rinsed, dried, weighed, and directly measured to obtain the coarse residue carbon pool. Soil samples below the sieve were added to a sodium polytungstate solution with a density of 1.8 g / cm³ and centrifuged at 1000 times gravity for 15 min. The floating matter was the free particulate organic carbon pool. The sediment was ultrasonically dispersed at an energy density of 300 J / mL and then subjected to repeated density separation. The second floating matter was the closed-state particulate organic carbon pool. The remaining heavy components were separated by 0.02 mm... The soil samples were separated into fine-grained heavy organic carbon libraries and coarse-grained heavy organic carbon libraries. Inorganic carbon was removed from carbonate-containing soil samples by hydrochloric acid vapor fumigation. The mass recovery rate of the five carbon libraries was greater than 95%, and the carbon recovery rate was 95% to 105%. The carbon 13 atomic fraction of each carbon library was obtained by dividing the difference between the carbon 13 atomic fraction of the labeled conditioning material and the carbon 13 atomic fraction of the same carbon library in the reference subregion by the difference between the carbon 13 atomic fraction of the labeled conditioning material and the carbon 13 atomic fraction of the same carbon library in the reference subregion. The labeled source carbon ratio was then calculated from the labeled source carbon ratio, the organic carbon content of the carbon library, the soil bulk density, and the soil thickness.

7. The method according to claim 6, characterized in that, Under no-till conditions, samples were taken on days 0, 7, 30, 90, 180, 360, 720, 1080, and 1440 after the addition of conditioning materials. Hydrothermal time was accumulated according to soil depth in the strips: data from 10-15 cm soil layer for shallow strips and data from 25-30 cm soil layer for deep strips. Hydrothermal time was the sum of the difference between the daily average temperature and 5°C for each day meeting the condition of a daily average temperature above 5°C and soil moisture content not less than 50% of field capacity. A model was established for the transformation of coarse residual carbon pools into free particulate organic carbon pools. The first-order kinetic equations, with rate terms representing the transformations of particulate organic carbon pools into closed-state particulate organic carbon pools, the transformations of fine-grained and coarse-grained organic carbon pools into fine-grained organic carbon pools, and the losses of each carbon pool in the form of mineralization and leaching, were used as the initial values, which were the measured values ​​on day 0. A nonlinear mixed-effect model was established using the natural logarithm of each rate constant as a parameter for fitting. Before fitting, the identifiability of the parameters was tested using the rank and condition number of the sensitivity matrix. The rank-deficient parameters were fixed one by one according to the Akaike information content criterion. The normalized trajectories of shallow and deep stripes were synthesized using calendar time as a common index and weighted according to their respective carbon input proportions; subsequent crops were sown on the day when the cumulative hydrothermal time in the 0cm to 20cm soil layer reached 240℃·d, and the latest sowing date was set.

8. The method according to claim 5, characterized in that, A full crossover experiment was conducted with three factors: the ratio of conditioning materials, the proportion of strips, and the tillage sequence. A total of 27 treatments were implemented, replicated three times, and started in three consecutive years, for a total of 81 plots. Alternatively, a D-optimal design was adopted, in which all interaction terms between the two factors could be estimated and there were no fewer than 24 different design points. The number of deep tillage operations in the three tillage sequences were all once in a four-year crop rotation cycle, but the year sequence was different. The soil depth for deep tillage was 30 cm, and the soil depth for rotary tillage was 15 cm. A blank subregion with the same proportion of strips was set up in each plot.

9. The method according to claim 8, characterized in that, Before the experiment, soil columns from 0cm to 60cm were collected from each plot, and the bulk density, cumulative soil mass, percentage of particles smaller than 0.02mm, organic carbon content of each carbon pool, and percentage of water-stable aggregates larger than 0.25mm were measured in 5cm layers as baseline values. Net carbon sequestration was calculated as follows: the minimum cumulative soil mass from 0cm to 45cm in all plots was used as the baseline soil mass, and the organic carbon storage was calculated using the equal-mass soil layer method. The increment of organic carbon storage in the treatment plot relative to the baseline value and the increment of organic carbon storage in the blank sub-region relative to the baseline value were taken. The difference in increments was converted to carbon dioxide equivalent at a ratio of 44 to 12, and then the carbon dioxide equivalent of the cumulative increase in methane emissions converted to 27, the carbon dioxide equivalent of the cumulative increase in nitrous oxide emissions converted to 273, and the emissions of the increase in diesel consumption for tillage operations converted to 2.68 kg of carbon dioxide per liter were subtracted. Greenhouse gas fluxes were measured every 7 to 10 days throughout the growing season, and were measured more frequently for 5 consecutive days after tillage operations, fertilization, daily precipitation greater than 10 mm, irrigation, and freeze-thaw events. Soil carbon dioxide fluxes were counted separately and used to verify mineralization and activation processes.

10. The method according to claim 9, characterized in that, The phase coordinate target values ​​are calibrated as follows: Net solid carbon is used as the dependent variable, and the independent variables are the linear and quadratic terms of mineral-bound recovery rate and closed-loop protection recovery rate, their interaction term, the main effects of conditioning material mass ratio and strip number ratio, their interaction term, carbon loading per unit strip length, and the random intercept of the start-up year. A response surface model is fitted, with the tillage sequence not included as a factor in the response surface model. Before fitting, a coordinate identifiability test is performed. The coordinate identifiability test is performed if the variance inflation factors of both mineral-bound recovery rate and closed-loop protection recovery rate are less than 2, the design matrix is ​​full rank and the condition number is less than 30, and the constraint optimum lies within the convex hull formed by the measured values ​​of mineral-bound recovery rate and closed-loop protection recovery rate. The total retention rate is the ratio of the sum of carbon reserves from marked sources in each carbon pool to the marked carbon input. The 95% confidence lower limit of the total retention rate is not lower than a pre-set value. Under the constraint of a certain threshold and the proportion of water-stable aggregates being higher than the baseline value, the constrained optimal point of the response surface model is obtained to obtain the target value of the phase coordinate. An alternative model with the tillage sequence as a factor is then fitted, and the alternative model is compared with the response surface model according to the Akaike information content criterion and cross-validation. When the mineral binding recovery rate and the closed-loop protection recovery rate do not simultaneously reach the target value of the phase coordinate within one crop rotation cycle, rotary tillage is adopted for all basic tillage operations within the same crop rotation cycle. The number of deep strips is selected based on the largest net carbon fixation amount predicted by the response surface model, under the constraint that the saturation deficit of fine-grained heavy component organic carbon is not less than 0. The saturation deficit is the difference between the estimated value of the fine-grained heavy component organic carbon capacity in the soil layer and the measured fine-grained heavy component organic carbon content, both expressed as the mass of carbon contained in each kilogram of soil. Finally, the conditioner composition, strip application parameters, target value of the phase coordinate, and tillage sequence are output.

Images