Modélisation du lessivage des bromures dans des cases lysimétriques par la méthode numérique des lignes Simulation of bromide leaching in pan lysimeters using the numerical method of lines

Résumé

La modélisation du transport des solutés dans un milieu non-saturé repose habituellement sur l'équation de dispersion-advection (EDA). Un modèle numérique (TSOL) a été développé en couplant l'EDA avec l'équation de Richards et en incluant le prélèvement de l'eau par les plantes. La résolution numérique a été effectuée par la méthode numérique des lignes (MNL) qui présente une grande simplicité de programmation et résulte en une très bonne précision numérique. La précision de TSOL a été testée avec les résultats d'un modèle d'éléments finis (HYDRUS), et avec des données expérimentales (profils de concentration et masses de bromures récupérés) collectées pendant 195 jours dans trois cases lysimétriques installées sur un sol non remanié cultivé en pommes de terre. La comparaison entre TSOL et HYDRUS montre que la solution de la MNL est similaire à celle des éléments finis. Toutefois, pour l'ensemble des cases et des profondeurs, les modèles ont montré une surestimation des valeurs de concentration avec un écart moyen entre les concentrations mesurées et simulées par TSOL variant de 22 à 112 mg/l. Pour les cases B et C, l'erreur moyenne de biais hebdomadaire entre TSOL et les masses de bromures récupérés, était d'environ 5 mg/semaine. Dans le cas de la case A, l'erreur moyenne de biais hebdomadaire de TSOL était de 39 mg/semaine.

Mots-clés:

• Modèle numérique,
• équation de transport des solutés,
• méthode numérique des lignes,
• case lysimétrique,
• lessivage

Abstract

Simulation of solute transport under transient unsaturated conditions is generally based on the dispersion-advection equation (DAE); a partial differential equation of the parabolic type under unsaturated conditions. The DAE has been solved by various numerical methods, such as finite elements and finite differences. However, these methods require advanced knowledge in mathematics and computer programming, in addition to specific adaptations to each problem in order to avoid numerical difficulties such as stability and convergence. The numerical method of lines (NML) can solve complex problems while keeping programming to a level accessible to a large number of engineers and scientists. The purpose of this article are (1) to develop and evaluate the numerical performance of a NML model (TSOL) that solves the DAE coupled with Richard's equation under unsaturated conditions; (2) to compare results of the TSOL model with those of a recognized finite elements model (HYDRUS); and (3) to validate the TSOL model with experimental data collected during 195 days under a potato field.

The experimental setup was installed on September 1994 in a potato field located at Saint-Pierre, Île d'Orléans, near Québec City. It consisted on three pan lysimeters (A, B, and C) with a surface area of 0.48 m² and a depth of 1.00 m, installed in an undisturbed sandy soil. On May 12, 1995, 15 g of KBr, dissolved in 60 ml of water, were applied uniformly over the surface of each pan lysimeter. The applied bromide was monitored until November 23, 1995. The monitoring period was divided into a first phase of 49 days, during which the soil was not cultivated and measures taken daily at 4 pm, and a second phase of 146 days during which measures were taken every Wednesday at 4 pm. The three pan lysimeters were sowded with potatoes on July 5 and harvested on September 5. Monitoring of the pan lysimeters included:

1. the drained water volume;

2. the water volume sampled by the pan lysimeter; and

3. the Br- concentration of all samples.

Numerical solution of the governing equations was obtained by the NML which belongs to the semi-discret methods consisting in discretising all independent variables except time, which is considered continuous for initial conditions problems. The discretisation of the DAE spatial variables was done by finite differences and resulted in a system of ordinary differential equations solved by LSODES (Livermore Solver for Ordinary Differential Equations Sparse); a solver used for systems with a sparse jacobian. The sparse nature of the jacobian results from our numerical procedure which solves simultaneously the DAE and Richard's equation.

The initial simulation time was fixed to May 12, 1995 at 4 pm and final time to November 23, 1995 at 4 pm with an hourly time step. The total depth of the pan lysimeter was simulated with a uniform internodal space of 1.0 cm. For each pan lysimeter, the initial pressure profile was measured by five pairs of tensiometers at depths of 7.5, 22.5, 45.0, 70.0 and 100.0 cm. The total mass of bromide applied was distributed equally over the first three upper nodes, and converted to concentration using the water content. At the soil surface, boundary condition for Richard's equation was taken as the hourly amount of rain fall (from planting to harvest), and as the net hourly water flux (from harvest until the end of the monitoring period). For the DAE, the boundary condition at the soil surface was of the third type. At the bottom of the pan lysimeters, a contant pressure head was assigned, which was the mean pressure measured by the deepest tensiometer; for concentration, a boundary condition of zero gradient was assigned.

For the two simulated variables (soil water bromide concentration and recovered mass of bromide), results of TSOL and HYDRUS were similar, showing an over-estimation of bromide concentration profiles but similar drained masses of bromide except lysimeter A, for which a large over-estimation was observed. This over-estimation may be explained by the presence of cracks between the soil and the plastic film surrounding the lysimeter. These cracks may allow a quick surface water flow along the sides during heavy rainfall. Because the solute was initially applied over all of the surface area of the lysimeters, a fraction of the solute might have migrated with the flowing water, and the remaining fraction by the soil matrix. The water flowing through cracks will quickly reach the bottom of the lysimeters. For all pan lysimetres and soil depths, the mean absolute error for weekly soil solution concentration profiles was 96 mg/l and the mean bias error 80 mg/l. For the mass of bromide recovered in lysimeters B and C, the weekly absolute mean error was 7 mg/week and the mean bias error 5 mg/week. For lysimeter A, the weekly absolute mean error and the mean bias error were the same, that is 39 mg/week.

For the simulation period, the numerical mass balance was negligible for the two models. However, the simulation time was longer for TSOL than for HYDRUS (18 versus 13 min). This difference is explained by the completely implicit differentiation scheme used by HYDRUS compared to the Backward Differentiation Formula used by LODES which is more complex, but frees the user from checking timestep precision. Considering the ease of programming and the resulting numerical precision, the NML has proven very effective in solving the solute transport equations in unsaturated conditions.

Keywords:

• Numerical Model,
• Solute Transport Equation,
• Numerical Method of Lines,
• Pan Lysimeter,
• Leaching