A computational model of the cereal leaves hydraulics

Environmental factors and plant architectonics significantly determine its water regime, namely, the water content and its movement through the tissues forced by the difference in water potentials, turgor pressure in cells, etc. In turn, the cumulative water regime affects the functioning and growth of cells, photosynthesis, and, as a result, the plant’s growth. The leaf contributes to the formation of the water regime and characterizes the contour of the plant’s adaptive system. Nevertheless, the data on the contribution of leaves to the total resistance to water transport in the plant and the structure of the hydraulic resistance of the leaf itself is still contradictory. This paper presents the formulation and justification of the monocots leaf hydraulics model based on Darcy’s law. The model was tested in computational experiments in the Comsol 4.3b package on idealized geometric models of leaf blades. Simulations showed the dependence of water potential distribution in xylem vessels and leaf mesophyll on the permeability of these tissues and on microclimatic parameters around the leaf. The adequacy of the model parameters selected as a result of testing is discussed.


Introduction
The plant water regime includes water content and its movement through plant tissues under the influence of a difference in water potentials, turgor pressure in cells, etc. It depends on environmental factors and plant architectonics and affects the functioning and growth of cells, photosynthesis, and ultimately, the plant's growth. It is known that plants grown in different climatic conditions differ in size, morphology, histology, and physiological parameters. The leaf makes a significant contribution to the formation of the water regime of the plant [1], and characterizes the adaptive system of the plant. Nevertheless, the contribution of leaves to the total resistance to water transport in the plant and the structure of the hydraulic resistance of the leaf itself remains controversial [2]. For example, the total hydraulic resistance of the leaf is summed up in equal proportions from the resistance of the leaf xylem and the resistance of paths outside the xylem in the mesophyll as emphasized in the studies [3,4,5]. In addition, it was argued that the dissipation of a significant difference in water potential at short distances within the leaf requires very high hydraulic resistance [6]. According to other authors, vascular resistance can range from 20% to 70% of the total leaf resistance [1].
Different classes of models are used for the theoretical study of the movement of water in a plant and its parts. Thus, the analogy between Ohm's law for electrical circuits and the linear dependence of the water flow rate on the pressure difference, which was observed experimentally using pressure chambers, led to the use of models with lumped parameters to describe the  [1] and links there). This type of model was also used to simulate the movement of water along a leaf [7]. Continuous models, or distributed parameter models, are another model type. Their theoretical basis is the models of convection and diffusion transport processes, including ones in porous media [8,9]. For example, [10] used models based on the Poiseuille and Darcy equations to study the movement of water in the xylem vessels and mesophyll of a plant leaf.
This paper presents a computational model simulating the hydrostatic potential distribution formation over the leaf of a monocotyledonous plant depending on the resistances of the xylem vessels and mesophyll. Preliminary results and their discussion on the example of a corn leaf were presented in the Comsol 4.3b package.

Water flows inside leaf tissues
Water flows inside the leaf tissues moving through the conductive system composed of a bundle of vessels, then leaves them and passes through the apoplast of the mesophyll, after which water evaporates into the air space inside the leaf, and in the vapor form exits through the stomata into the atmosphere. The driving force of the water flow is the gradient of the hydrostatic potential [11,12].
In the mesophyll, the emitted from the xylem water fluxes can follow along three pathways. Firstly, water unavoidably moves along the apoplast formed by the interconnected space of the cell walls and finally exits into the intercellular air space inside the leaf. In addition, water can pass through the cells' cytoplasm, entering and exiting through the pores and plasmodesmata (second path) and the aquaporin water channels (third path) [13,14]. The throughput of the second and third paths is known to depend strongly on temperature and other external factors [13].

Model equations
While modeling the movement of water inside the leaf, both conducting vessels and mesophyll tissues can be considered as porous media with different geometries and scales of the liquid phase moving in a space that is free from the solid phase of the medium (for example, from macromolecular complexes, etc.). This assumption allows us to apply the Darcy model to all leaf tissues [15].
The model is formulated as a following system of two partial differential equations. Darcy's law:ū where u is the Darcy velocity (SI unit: m s −1 ), κ is the permeability of the porous medium (SI unit: m 2 ), µ is the fluid's dynamic viscosity (SI unit: P a s), p is the fluid's pressure (SI unit: Pa), ρ is fluid's density (SI unit: kg m −3 ), g is the magnitude of gravitational acceleration (SI unit: m s −2 ), ∇D is a unit vector in the direction over which the gravity acts. Here the permeability, κ, represents the resistance to flow over a representative volume consisting of many solid grains and pores. When modeling a leaf, the contribution of gravity to the water potential is assumed to be zero. Continuity equation: which in the case of incompressibility of the liquid takes the form: where ρ is the fluid density (SI unit: kg m −3 ), ϵ is the porosity defined as the fraction of the volume occupied by the pores, Q m is a mass source term (SI unit: kg m −3 s −1 ). Substituting Darcy's law into the continuity equation, we obtained a governing partial differential equation for the pressure function p(t, x, y, z):

Model testing
Adequacy testing for the model consisted of evaluating and optimizing the model's parameters so that the solution to the model satisfies some optimization criteria consistent with the experimental data. We used transpiration rate, the average hydrostatic potential of water in the leaf, and the contribution of xylem to total leaf resistance as optimization criteria.
Numerical solutions for the model of monocotyledonous leaf hydraulics were obtained in the Comsol 4.3b package (Information and Computing Center of Novosibirsk State University).

Geometric model of the leaf
The data on the size and structure of the vascular system of the leaf were obtained on a sample of a maize plant ( textit Zea mays L.) cultivar B73. The model leaf had the following characteristics: length 80 mm, maximum width 10 mm, thickness 0.19 mm, the diameter of large vessels 70 µm.
We built a simplified three-dimensional geometric model of the leaf representing an elongated ellipse with a truncated base and xylem vessels located along the longitudinal axis. Since the leaf length and width values are 2-3 orders of magnitude greater than the thickness, in a rough approximation, the distribution of variables and parameters over the thickness is homogeneous. The main simplifications consisted of the following points: (a) only the conductive bundles of the first-order xylem vessels were considered (the smaller vessels are taken into account in the mesophyll parameters), (b) the conductive bundles of the xylem vessels in the cross-section were modeled as rectangular regions with dimensions of 6 µm and 190 µm (leaf thickness). Simplified models of the leaf and its fragments were constructed using the geometric modeling tools in the Comsol 4.3b package.

Stationary problem formulation
To study the parameters of water movement along xylem vessels and the mesophyll inside the leaf, we formulated a stationary problem with boundary conditions. Since when measuring the parameters for the leaf water regime in short-time periods, we can neglect the absorption and release of water in the processes of respiration and photosynthesis and assume that Q m = 0, then the equation 4 takes the following form: We took the following boundary conditions (BC): BC 1: There is no flow at the lateral border of the leaf area and at the ends of the vessels near the leaf apex −n · ρū = 0, wheren is normal vector to the region boundary. BC 2: Pressure at the vessels borders in the leaf base is p σ i = −4 · 10 3 P a, where σ i ∈ Σ(x = 8) are border areas belonging to xylem vascular bundles. This pressure is required to raise the water column to a height of 40 cm. We assume that the leaf is at this height. The water flow through the lower and upper surfaces of the leaf is equal to the flow of the water vapor through the stomata, and is described by the formula: where g v (SI unit: kg m −2 P a −1 s −1 ) is the conductivity of the vapor flow from the leaf into the atmosphere through the stomata and then through the boundary layer near the leaf surface; a = 0.611 kP a, b = 17.502, c = 240.97 • C are empirical constants suitable for ecological-physiological tasks; θ is the Celsius temperature; h a is relative humidity, M w = 1.018 kg mol −1 is molar mass, p (SI unit: P a) local water pressure at the leaf border, The equation 6 is obtained from the following reasoning. The dependence of the saturation pressure of water vapor over the free surface of distilled water on the surface temperature satisfies the following relation: where θ is the Celsius temperature, a = 0.611 kP a, b = 17.502, c = 240.97 • C are empirical constants suitable for ecological-physiological studies [12,16]. If the water is not free and/or not distilled (for example, capillaries or solution), then the relative humidity of water vapor above the surface is described by the formula: where h r is relative humidity, M w = 1.018 kg mol −1 is molar mass, ψ (SI unit: J kg −1 ) is water potential, R = 8.31 J K −1 mol −1 is universal gas constant, T is the Kelvin temperature [12]. Since in the Comsol package the water potential is included in Darcy's law in the form of pressure, we rewrite the equation 8 in the following form: where ρ (SI unit: kg m −3 ) is density of water. Within the linear model, the vapor flow from the leaf to the atmosphere is determined by the following relationship: where e a (SI unit: P a) is vapor pressure in the atmosphere e L (SI unit: P a) is saturation pressure in the leaf, g v (SI unit: kg m −2 P a −1 s −1 ) is conductivity for vapor flow from the leaf through the stomata and then through the boundary layer near the leaf surface into the atmosphere. The vapor pressure in the atmosphere and the leaf temperature (here it is taken equal to the outside air temperature) are external parameters for the model. Local vapor saturation pressures in a leaf at a given temperature and local water potential are calculated using the formulas 7 and 9; substituting 7 and 9 into 10, we get the formula 6.

Using a continuous model for estimating the parameters of a lumped model
Depending on the modeling goals while simulating the water movement through a plant and its parts, both a model with distributed and lumped parameters are used. That is why it is important to associate the parameters of these model types. In particular, lumped-parameter models make it possible to plan and interpret the experimental studies of the plant water regime. The correspondence between the concepts and parameters of the Darcy model and the model similar to Ohm's law are discussed in [17]. In this work, the ratio of the xylem resistance to the total leaf resistance for water flow was one of the criteria for the model parameters optimization.
The resistances of xylem and mesophyll were calculated using the following formulas. Within the framework of linear models, the flow Φ (SI unit: kg m −2 s −1 ) through the domain and the generalized force ∆Ψ (N ) on the boundary of the domain are connected through the ratio: where R is flow resistance. According to the conservation law, the flows of water entering the region of vascular xylem bundles and leaving the mesophyll region are the same and equal to the total evaporation. The flow of water through each area can be calculated, for example, using the formula: where Σ l is the boundary of the mesophyll region, through which the water leaves. For conducting xylem vascular bundles ∆Ψ is the difference between the integrated pressures along the parts of the boundary where the flow exits and the pressures at the entrance of the vessels: where Σ x i, Σ x o are the boundaries of the vessels where the flow enters and exits, respectively. For the mesophyll area, ∆Ψ can be calculated by the formula: where p(σ) is the pressure distribution along the border of the leaf lamina. From the formulas (13)- (14) forses were calculated to be equal to ∆Ψ x = 6 N and ∆Ψ m = 100 N , respectively. This leads to the fact that the contributions of the xylem vessels and mesophyll resistances to the total leaf water flow resistance are estimated to be approximately 6% and 94%, respectively. Many studies note that experimental data on the contribution of xylem to the total resistance of the leaf to water flow are still contradictory. For example, the vascular resistance can be estimated to be from 20% to 70% of the leaf resistance [1]. Note that the results of model calculations based on experimental data lead to an estimate of 2% [3,18]. Such an underestimated result for xylem resistance maybe since the resistance of the bundle sheath cells was not taken into account.

Estimation of the model parameters and their optimization
The estimate for the permeability coefficient for a bundle of xylem vessels in the Darcy model was made based on an analogy between Darcy and Poiseuille laws. From the relationsū = κ µ∇ p andū = d 2 32 µ∇ p we conclude that κ = d 2 32 [20]. According to this formula, with a diameter of one vessel equal to 5.66 µm the permeability coefficient for a bundle of xylem vessels is κ x = 10 −12 m 2 .
The estimation of the permeability coefficient for the apoplast material was obtained from the experimental data on the specific water conductivity [21], L pp = (1.0 ± 0.6) 10 −3 cm 2 s −1 bar −1 . From the definition of the flow through the cell wall, we can conclude that q = L p (∆P ), where L p is water conductivity [11]. According to Darcy's law, the flow through this wall is q = κ µ 1 δs (∆P ), where δs is cell wall thickness. Considering that L p = Lpp δs , we can calculate that κ m = L pp µ ≈ 10 −15 m 2 .
When calculating the movement of water in the mesophyll along the apoplast, it is necessary to consider the ratio of the thickness of the cell wall and its average diameter. Depending on this, the coefficient of permeability of the region in the model corresponding to the mesophyll can be two or three orders of magnitude less than the permeability of the material. On the other hand, the effective permeability of the region may be greater than the above value if the water fluxes occur both along with the mesophyll cells and along with the air space inside the leaf, which can be partially filled with water, depending on the moisture content. These variations may reflect an ecophysiological reality.
To optimize the model to the experimental data, we varied the permeability coefficients for xylem and mesophyll and the conductivity of the stomata. Criteria for selecting suitable values were published estimates of the evaporation rate, average water pressure in leaf tissues, and the contribution of xylem and mesophyll to the resistance of the leaf to water flow. At the same time, we tested the sensitivity of these indicators to parameters variations. As known during the day, the loss of water for transpiration can be several times higher than the weight of the plant. The exact coefficient depends on many external and internal factors and estimates the transpiration rate. We optimized the values of xylem and mesophyll permeabilities at a transpiration rate 5 · 10 −6 kg m −2 s −1 , which corresponds to 4 times the mass of the leaf. Figure 1 shows the dependences of the xylem and mesophyll resistances on the permeability of the mesophyll. It is interesting to note the nonlinear dependence of the resistance of the mesophyll on its permeability (figure 1 B). The reason for this may follow from the distribution of the hydrostatic potential in the leaf. of the conductivity variation interval. The change in the pattern of the isolines allows us to conclude that the water flow lines in the leaf are changing (the flow vectors are perpendicular to the isolines). Taken together, this leads to a nonlinear dependence of the resistance of the mesophyll on its permeability.
Xylem permeability affects both xylem resistance and mesophyll resistance as seen in figure 2. The reason for the effect of xylem permeability on the resistance of the mesophyll is a change in the pressure distribution in the mesophyll with a change in the pressure in the xylem (figure 2D, E).
As mentioned above, we choose the following optimization criteria: (i) the transpiration rate should be such that the mass of evaporating water per day exceeded four times the mass of the leaf; (ii) the average hydrostatic potential of water in the leaf is about 5 · 10 5 P a, and (iii) the contribution of xylem to the total leaf resistance exceeds 10%.
The values of mesophyll and xylem permeability κ m = 10 −17 m 2 and κ x = 10 −4 m 2 , respectively, satisfy this set of criteria. It turns out that the selected xylem permeability differs by two orders of magnitude from the initial estimate. Similar discrepancies are noted in [9], where, as possible reasons, the authors note the insufficient knowledge of the fundamental issues concerning the water regime of plants.

Influence of the transpiration rate on the leaf water regime characteristics
The equation 6 shows the dependence of the transpiration rate on stomatal conductance, air and leaf tissue temperature and on the difference in hydrostatic vapor potentials within the leaf and the atmosphere. The conductivity of stomata depends on the morphological parameters and the physiological processes in the plant. The difference in hydrostatic potential depends on  Figure 3. Influence of the transpiration rate on the average and minimum values of the hydrostatic potential.
the temperature and water potential of the liquid phase of the water in the leaf. However, at absolute values of the hydrostatic potential of the liquid phase from 4 · 10 4 to 3 · 10 6 P a, the hydrostatic potential of water vapor inside the leaf differs from zero by a few tenths of a percent, which makes it possible to take the relative humidity of the air inside the leaf equal to one.
On the highly simplified level, varying the transpiration rate could represent the aggregate effect of all external and internal factors on the leaf water regime. Figure 3 shows the results of varying the transpiration rate in the range from 2.4 · 10 −6 to 9.8 · 10 −6 kg m −2 s −1 , which corresponds to a change in the relative air humidity from 20% to 80%.
Note that change in the hydrostatic potential in the extracellular space, which is not taken into account in the current version of the presented model, determines the flow of water in the leaf and also directly affects the exchange of water in the cell and therefore affects the physiological processes in the leaf.

Conclusion and perspectives
Common biophysical concepts of plant cells and tissues and data on the ecophysiological processes of water flow in the plant-atmosphere continuum formed the basis for a threedimensional model of water movement through a leaf into the atmosphere. The simulation results conclude that the model considers the determining processes and leads to adequate parameter estimates. The model makes it possible to obtain plausible values for the distribution of water potentials (pressures) and to estimate the contributions of xylem and mesophyll to the leaf water flow resistance depending on the essential parameters of the plant microclimate.
The model will allow planning experiments to study the water regime of a plant leaf, including MRI, and interpreting their results. In the future, the model can also include calculating the leaf temperature (thermal conductivity and various types of mass transfer), which is established as a result of radiation exchange and other forms of energy exchange. The model can serve as a component for a more general model of the water regime of the plant.
At the same time, computational experiments have shown the relevance of developing a series of different scale models of the leaf water regime for studying the contribution of leaf cell architecture to the parameters of the large-scale model presented in this work.