Multiscale Modeling of Pollutant Uptake by Mangroves

Models for pollutant uptake by mangrove plants are developed at different scales and applied to the environmental fate of a pollutant in an estuary. Employing cohesion-tension theory, a 3-dimensional model of water and substance flow in young mangrove trees is set up in form of porous media equations. Water transport is conceived as a continuous hydraulic process driven by canopy transpiration. State variables are water potential and pollutant concentrations in the soil, roots, xylem, core and canopy. At catchment scale, the mangrove forest is conceived as a flow reactor. Upscaled models are derived from a 3-dimensional single plant model by fitting a compartment model in form of ordinary differential equations to data obtained by spatial integration over the domains of the 3-dimensional model. These equations are imbedded as reaction terms into the shallow water equations for riverine transport. The model is applied to the dispersal of pollutants in an estuary.


INTRODUCTION
Industrial pollutants such as heavy metals and organic substances discharged into coastal areas constitute a major threat to the environment and ecosystem.For tropical and subtropical regions, mangrove ecosystems are dominant in the intertidal zone.Their ecological and socioeconomic values have been widely recognized.Mangroves also retain pollutants from wastewater and they can serve as a means for immobilization and removal of pollutants [1,2].
The ability of mangroves to retain pollutants can be used in applications in a variety of technologies such as waste water treatment, biotechnology, and especially in emergency response in cleaning up polluted rivers on the catchment level by constructed wetlands [3,4,5,6].For the assessment of such methods the determination of the capacity and threshold of mangrove trees to uptake pollutants need to be specified first at the level of single trees by analysing the underlying plant physiological mechanisms.In a second step, these processes have to be related to the landscape level.This concerns the flow dynamics of pollutants in a catchment and their interaction with mangrove trees with respect to uptake and storage and toxic effects on plant growth.The motivation of this work is derived primarily from the need to bridge the micro (single trees) and macro (catchment) scales in the sense that one might still be able to capture the macroscale behaviour of the system with the help of the microscale models.
Any system of interest can always be described by a hierarchy of models of different complexity i.e., both physical and biochemical processes take place in more than one scale level and the description of each level is determined by the processes that occur in the lower ones.By the aggregation of models at the micro scale and integrating them into models at the macro scale, we may be able to describe the hierarchical system in a consistent way without too much loss of accuracy [7,8,9,10].
We present models for simulating the propagation and storage of pollutants in mangrove trees at different scales.This concerns both the application of multiphysics in biology and the upscaling approach from single trees to landscape scale.The uptake and transport of pollutants in single trees constitutes a typical multiphysics problem involving physical and biological processes.The transport of water and matter in the soil and in the plant can be described in the frame of a physical theory based on the cohesion-tension concept, whereas the control of evaporation by the canopy is a biological process involving complex biochemical mechanisms.
The second aspect concerns the modelling of the environmental fate of pollutants at larger scales e. g. catchment scale.This involves the upscaling of single tree models to a mathematical form which can be coupled to large scale dispersal models.Here, we show how 3-dimensional models for single trees in form of partial differential equations can be aggregated to systems of ordinary differential equations.The aggregated model can then be implemented into dispersal models, here the shallow water equations, as reaction terms.The models should ultimately support the design of phytoremediation systems and the landscape planning under the aspects of phytoremediation.The research approach is shown in Figure 1.
This paper focuses on the aggregation of small-scale models and implementation into large scale transport models enables the simulation of the propagation and retention at landscape scale.

THREE DIMENSIONAL MODEL FOR WATER AND POLLUTANT TRANSPORT IN YOUNG MANGROVE TREES
The governing equations are based on the cohesion tension theory.The concept of this theory is based on the water potential.The water potential gradient between the soil and the leaves is the driving force for water transport (cf. Figure 2).Water transport in the soil and tree is conceived as a continuous hydraulic process, which is modeled by porous media transport equations [12,13,14,15].State variables are water potential and contaminant concentrations in the soil, roots, xylem, core and canopy.For each domain, a porous media equation (Richards equation) is set up which is coupled to adjacent domains via internal boundary conditions.
[ ] Contaminant transport equations are coupled to the water transport equations via the Darcy velocity and the dispersion tensor.
In all domains, contaminant species can be adsorbed or bound to soil particles or tissue respectively.
The biological sub model [17] for the transpiration is presented in Box 1 for a better understanding of the interrelations between the environmental variables and the state variables of the leaf.The transpiration rate is driven by the difference in specific humidity at the leaf surface and the saturation specific humidity of the air and the wind speed.The bulk transfer of latent heat is controlled by the stomatal conductance, which in turn depends on photosynthetic radiation, leaf temperature, vapor pressure deficit and leaf water potential.

Geometry
The geometry of the plants were derived from the young mangrove trees used in our experiments.Mangrove trees were planted into containers with holes at the side and at the bottom and placed into a water basin, which was flooded twice a day in order to simulate the tide conditions.These experimental conditions are implied in the boundary conditions below.The geometry model was created in MATLAB and implemented into the software COMSOL MULTIPHYSICS.
Box 1: Biological submodel for the control of transpiration.For further explanation, see [15]. Figure 3: Geometry and meshing for young mangrove trees in containers with the domains soil, stem, branches and leaves.To capture the variability of the plants, several geometries were generated differing with respect to leaf numbers and tree size.

Initial and boundary conditions
The flow between different domains is mediated by the water potential difference and by the concentration difference respectively.Active transport is not considered here.
) ( External boundary conditions concern the infiltration of water and pollutants into the soil and the water loss by evaporation.We consider time varying boundary conditions as motivated by the experimental setup.Containers with young mangrove plants are placed into a basin which is flooded twice a day simulating the conditions in a tidal zone.Water infiltrates the containers via the surface and by holes in the container.

Some simulation results
The boundary value problem was implemented into the software COMSOL MULTIPHYSICS [18] using the coefficient form of the PDE interface.COMSOL is a finite element analysis solver and simulation software, Figure 4 shows snapshots of the water flow and concentration distribution in the soil-plant system.The snapshot of the water flow in Figure 4a was taken in the early phase of the inundation, when the container is not yet totally immersed.So one clearly distinguishes the water flow towards the surface due to evaporation and the water flow to the plant root and subsequently to the leaves driven by transpiration. Figure 4b shows the distribution of the pollutant in the soil and in the plant after 4.5 hours.
More details of this single tree model and simulations can be found in [15].

UPSCALING
The basic strategy for implementing an upscale procedure is to link together simulations of the lower level models on small spatio-temporal domains in order to mimic the behaviour of a system at large scale.Assume that the form of the large scale model is known, but some of the details of the model are missing.In this case, an aggregated version of the microscale model can be used to supply whatever data are needed in the macro-model.For the simulation of the environmental fate of pollutants at landscape scale, pollutant uptake and transport in single plants have to be upscaled.Therefore on the landscape scale the detailed 3-dimensional transport model for individual trees was replaced by a compartment model in form of ordinary differential equations.

Outline of the procedure
The following steps are involved in the upscaling procedure

Construction of compartment models
In a first step, the 3-dimensional model in form of partial differential equations is replaced by an aggregated model in form of ordinary differential equations.To this end, a compartment model is devised where each compartment represents a spatial domain of the 3-dimensional model (cf. Figure 5).This implies that concentrations are homogeneous in each compartment, which is regarded as a "well stirred" medium.So the spatial information within a compartment is lost.The only spatial information remaining rests upon the allocation of a variable to a compartment.In each domain, the substance is partitioned between the liquid and solid phase, so there are six state variables for the aggregated model.These are either the concentration or the mass of the substance.To construct aggregate models, numerical experiments with the 3dimensional model are performed.Initial and boundary conditions have to be kept simple to facilitate parameter and model identification.For each domain, the total amount of the substance is evaluated by integration of the concentration field over the respective domain.

∫∫∫
The resultant time series "data" are then used for model identification and parameter estimation.
2. Fit of these "data" to a compartment model It is obvious that aggregated models cannot match precisely the data generated from the 3dimensional model.There is always some bias.However, these data are completely deterministic.In nature, trees even from the same age class are different in size and physiological state.So in principle, one could generate data with a stochastic component by performing numerical experiments for a population of 3-dimensional models with different geometries (cf. Figure 3).Here, we employ a simpler procedure.The time series data generated by the model of one tree are superimposed by normally distributed errors.The standard deviations are chosen under the assumption of a constant coefficient of variation.Parameters are then estimated by nonlinear regression procedures for systems of nonlinear ordinary differential equations.This is achieved by embedding a solver for initial value problems (ODE 45 solver in MATLAB based on an explicit Runge-Kutta (4,5) formula) into the nonlinear optimization algorithm fminsearch in MATLAB, which uses the simplex search method of Lagarias et al. [19].We define an aggregated model as appropriate if the curves are matching the generated data within certain limits (two standard deviations).Figure 6 shows an example.
In the next step, spatial information at the higher scale (landscape, catchment) is reintroduced into the compartment model by spatialization of the model parameters.For example, the kinetic parameter determining the rate of uptake from soil to roots depends on the surface of the root which is related to the biomass and thus to the age of a tree and on the density of the trees.In this way single trees are replaced by a spatially distributed "reactor".

Combining the reactor model with a large scale transport model
In the last step, the spatialized reactor model is coupled to a water transport model for rivers and catchments.Technically spoken, the biological part of the model now consists of the reaction part in a system of partial differential equations.

Compartment model
In the aggregated model the state variables are the mass of the substance in the domains, in the spatialized model their dimension is density (mass/unit area).
In each domain, the substance is partitioned between the liquid and solid phase, so there are six state variables for the aggregated models.A compartment scheme as presented in Figure 5 was devised and modelled by the following system of linear ordinary differential equations representing the mass balance.
Equations (15 and 17) are the mass balance equations for the liquid and solid phase in the soil respectively  with input function F(t) Equations (20 and 21) are the mass balance equations for the liquid phase and tissue in the leaves respectively.To assess this model and to estimate parameters, time series of the amounts of the substance in the domains soil, stem including roots and leaves of the 3-dimensional model were generated.This was done for a simple batch experiment with a constant infiltration rate at the upper boundary, drainage outflow at the lower boundary and a mean value of the photosynthetic active radiation.Figure 6 shows the resultant fit to the deterministic data with superimposed normal errors.The aggregated model matches well the time series of the integrated concentrations generated by the 3-dimensional model.Although the partial differential equations of the 3-dimensional model are nonlinear, the linear compartment model matches the data quite well.However, this was proven only for the simple initial and boundary conditions.

Embedding into a real landscape
For upscaling, the model parameters have to be spatialized with respect to the stand density of the mangrove forest and to the biomass of the trees.Since in the Can Gio reservation the mangrove trees were reforested, the forest is made up of areas with the same age cohorts.The dependence of the model parameters on the biomass which is related to age are derived from the following growth model for mangrove trees, which was parameterized by Nguyen [20].The state variable of the growth model, y, is the diameter at breast height (dbh), which is usually employed in forestry.The factors fi(ei) are normalized to one and model the influence of environmental variables on the growth rate such as salinity or stand density.Figure 7 shows growth curves under several combinations of the environmental factors "salinity" and "stand density".The two dimensional shallow water equations are coupled with the convection diffusion equation for pollutant transport and the reaction equations [21].The shallow water equations constitute the conservation of mass (Equation 26) and the conservation of momentum (Equations 27 and 28).
The depth averaged pollutant transport is coupled with the shallow water equations ) , The pollutant concentration is coupled to the above aggregated model with the spatialized parameters (Equations 32 -37).
This model is specified for a catchment via the bathymetry f(x,y) and Manning's roughness coefficient n, which reflects the distribution of the vegetation.

Initial and boundary conditions
At the upper boundary, the height z was kept fixed, at the lower boundary z was varied according to the tide.The pollutant was discharged from a point source either continuously or as batch input.

APPLICATION TO A SECTION OF THE THI VAI RIVER IN VIETNAM
Thi Vai River is situated in the southern part of Vietnam separating an industrial area at the eastern bank from the Can Gio mangrove forest reservation in the west.One of the major pollutants in this region is chromium.Mangroves are highly efficient absorber of pollutants.
As our experiments [22] and field studies have shown mangroves have a high potential for phytoextraction of metals such as chromium.The 2-dimensional reaction and transport model was applied to the dispersal and storage of chromium in a section of the Thi Vai river.For this purpose, a simplified analytical bathymetry model was used.Two scenarios were considered: i) direct discharge from a point source into the river ii) discharge into a mangrove wetland prior to discharge into the river.
Figures 8a-c show snapshots of the dispersal of the pollutant in the case of direct continuous release into the river.-c show the respective snapshots for the case that the pollutant is released into an area densely covered with mangroves serving as a waste water treatment plant prior to the release into the river.For the latter scenario concentrations in the river are reduced and the mangrove forests at the river bank are shielded from pollution (Figure 10).Water flow in the enclosure with the mangroves is retarded due to the high resistance in the planted area, therefore the residence time of the pollutant is increased leading to higher concentrations in the soil (cf. Figure 12) and subsequently in the mangroves.Figure 12 shows the time course of concentrations of the pollutant in mangrove plants at a point in the center of the enclosure.On the basis of these simulations we suggest the construction of a pilot wastewater treatment plant on the basis of mangrove trees.

DISCUSSION AND CONCLUSION
Multiphysics modelling is not limited to physics.We have demonstrated that Multiphysics can be applied to biological systems combining physical and biological subsystems.However, there are differences in the governing equations for both realms.Whereas in physics models can be derived from first principles such as the conservation of mass or momentum, in biology frequently heuristic approaches prevail.This can be seen in the evaporation sub-model (see Box 1).The stomatal control function is composed of empirical functions describing the influence of temperature and other environmental variables on stomatal opening and thus on the rate of transpiration.The model equations both for the physical and for the biological part are highly nonlinear.Because of these nonlinearities and the extremely high differences in water potential between air and leaves, boundary and initial conditions and net size and structure must be set up very carefully in order to achieve convergence.The experience with this model shows that it is feasible to add more biology, e.g.photosynthesis, transport of assimilates in the phloem and, in the case of mangrove trees, the extrusion of salt.
Another important issue is the upscaling of complex models in form of partial differential equations.We have devised a simple procedure for the derivation of aggregated models in form of ordinary differential equations: treating the complex model as a real system, performing numerical experiments and identifying compartment models.It should be mentioned that only linear compartment models were identified although the underlying 3dimensional model includes nonlinearities.Also, the model still pertains to a single tree.
In order to proceed to large scale models, the parameters of the local aggregated model have to be spatially distributed and must correspond to the local vegetation cover.We have demonstrated how this can be achieved for a mangrove forest.As a practical application, the effect of a phytoremediation system on the dispersal of a pollutant in an estuary was simulated.The simulations have shown that the building up of a phytoremediation system can reduce the pollution in the downstream area.Further developments should imply parameter estimation in models in form of systems of partial differential equations, addition of more biology to the 3-dimensional model, possibly the derivation of nonlinear aggregated models and the coupling of the landscape model to a geographical information system.

Figure 1 :
Figure 1: Research concept: Field data and the data from the operation of a pilot plant for phytoextraction of metals are the basis of the model developments.

Figure 2 :
Figure 2: Cohesion tension concept.The water potential gradient is the driving force for water transport in the soil-plant system.Figure adapted from Kirkham[11]

Figure 3 :
Figure 3: Snapshots of water flow (a) and distribution of a pollutant (b) in the soil plant system.

Figure 4 :
Figure 4: Compartment scheme for matter transport in a tree.The parameters kij are kinetic constants.
19) are the mass balance equations for the xylem and core respectively.

Figure 6 :
Figure 6: Fit of the aggregated model to data from the 3-dimensional model.The data are superimposed by normal distributed errors.

Figure 5 : 4
Figure 5: Growth curves of mangrove trees under different combinations of the environmental factors "salinity" and "stand density ".

Figure 6 :
Figure 6: Snapshots of the dispersal of a pollutant in a section of the Thi Vai river estuary.The pollutant is directly discharged into the river.

Figure 7 :
Figure 7: Snapshots of the dispersal of pollutant which is released into an enclosure which is densely covered with mangrove trees prior to the discharge into the river.

Figures 9a
Figures9a-cshow the respective snapshots for the case that the pollutant is released into an area densely covered with mangroves serving as a waste water treatment plant prior to the release into the river.For the latter scenario concentrations in the river are reduced and the mangrove forests at the river bank are shielded from pollution (Figure10).

Figure 8 :Figure 9 :
Figure 8: Distribution of the pollutant in mangrove trees if a) the pollutant is directly released into the river b) the pollutant is released into the enclosure with mangroves prior to the discharge into the river

Figure 10 :
Figure 10: Time course of concentrations in mangrove trees in the center of the enclosure.