Interactive comment on “ Technical Note : Simple formulations and solutions of the dual-phase diffusive transport for biogeochemical modeling

This paper introduces an efficient scheme for dual-phase diffusive transport simulation. This approach has the advantage of being easily integrated in large-scale distributed environmental models. The paper is well written, and I recommend publication as is. The case studies are carefully designed, the math is well developed and clear, and the numerical simulations and benchmark experiments against analytic solutions are convincing. Of course, one can request additional benchmark tests involving transport in two or three-dimensions yet this is beyond the scope and method of solution of most climate models.


Introduction
The interest in predicting fluxes of various biogenic greenhouse gases and their interactions with climate change has motivated the development of many terrestrial biogeochemical models; e.g., ecosystem methane models (Walter and Heiman, 2000;Zhuang et al., 2004;Tang et al., 2010;Riley et al., 2011), nitrification-denitrification models (Venterea and Rolston, 2000;Maggi et al., 2008), water-CO 2 isotope mod-els (Riley et al., 2002), and generic reactive transport models that attempt to integrate as many biogeochemical processes and chemical species as possible (e.g., Simunek and Suarez, 1993;Grant, 2001;Tang et al., 2013).To resolve the depthdependent dynamics, these models in general represent multiphase (aqueous and gaseous phase) diffusion processes and often assume negligible advection.
The equation for multiphase diffusion has been represented in different forms by different authors (Table 1).However, the numerical implementation of the equation is often vaguely described (either by referring to other publications or by mentioning the numerical scheme) or is convolved with other technical details, making the model difficult to understand or replicate by other researchers.In many cases, however, one does not need to represent all the processes typically included in a complicated reactive transport model to understand a particular problem.For instance, when soil moisture and temperature data are available together with soil respiration, one only needs a diffusion model to evaluate belowground CO 2 dynamics (Davidson et al., 2006).Therefore, ecosystem models would benefit from a simple mechanistic formulation and numerical implementation of the dualphase diffusion problem.
In this note, we categorize the existing formulations of the dual-phase diffusion problem into three forms, and recommend one that can be most easily implemented numerically in land models.We hope that this effort will help researchers who wish to develop simple but mechanistic transport models for their particular problem.
Published by Copernicus Publications on behalf of the European Geosciences Union.Walter and Heimann (2000), Zhuang et al. (2004).
The gaseous phase is used as the primary variable.The model assumes gas diffusion to dominate in unsaturated soil and aqueous diffusion to dominate in saturated soil.The water table is assumed to be at the layer interface.Diffusivity varies continuously with soil moisture.Venterea and Rolston (2000), Riley et al. (2011) ∂C ∂t = ∂ ∂z D ∂C ∂z + S Bulk tracer concentration is used as the primary variable.Diffusivity varies continuously with soil moisture.Special care is put to the water-air interface.Tang et al. (2010Tang et al. ( , 2013)).
Tracks gaseous and aqueous phases of a given tracer separately.Gas dissolution and exsolution are considered explicitly.Diffusivity varies continuously with moisture.Maggi et al. (2008) 2 Methods

Governing equations
In this section we derive the relevant mass balance differential equations from first principles.Throughout this note we assume advection is treated using the operator splitting method (Tang et al., 2013), or is negligible.Considering the diffusive mass transport problem (Fig. 1), the dual-phase diffusive flux from layer j − 1 into j is where subscript w indicates aqueous diffusion and subscript g indicates gaseous diffusion.Soil moisture is represented by θ and air-filled porosity by ε.The fluxes (F j −1→j ) are imposed at the upper and lower boundaries of layer j .Tracer concentrations are designated by C with appropriate subscripts.In defining the effective aqueous (D w ) and gaseous (D g ) diffusivities, we assume the tortuosity has been considered appropriately (e.g., Moldrup et al., 2003).A full list of symbols is given in Appendix Table A1.
Applying the law of mass balance to layer j and in the limit of small z j and t, one obtains Figure 1.Schematic of the dual-phase diffusive transport problem: solid lines are the interfaces between control volumes and the dashed line is the center of the control volume.All symbols are defined in Appendix Table A1.
where S (C, z) is the tracer source due to processes other than diffusion and C is the bulk tracer concentration including both gaseous and aqueous phases.From the assumption of equilibrium between aqueous and gaseous phases (e.g., Tang et al., 2010): where α is the Bunsen solubility coefficient, one obtains, from substitution of Eq. ( 3) into the diffusion and temporal derivative of Eq. ( 2b): By further defining a bulk diffusivity and assuming one finds for the bulk tracer: Clearly, Eq. (4a) (gas-primary form), Eq. (4b) (aqueousprimary form), and Eq. ( 7) (bulk tracer form) are equivalent, but Eq. ( 4) are more convenient to solve because the tracer sources are in general given as aqueous reactions or gaseous sinks (e.g., plant transport, ebullition) and the necessary phase conversion can be done easily through Eq. ( 3).In particular, Eq. (4a) (the gas-primary form) is the most convenient for simulating volatile tracers and can be applied to variably saturated soil without the need for special care of the air-water interface, as was done in a few existing wetland-CH 4 models (see remark in Table 1).We note that Eq. ( 4a) was also used by Simunek and Suarez (1993) to model CO 2 transport in soil.
At the top boundary, conditions are in general given as gas tracer concentrations, which are connected to the gas concentration of the top numerical layer as where r a is atmospheric resistance and r s is soil resistance (see Tang and Riley, 2013 for a detailed discussion).At the bottom boundary, zero flux conditions are often imposed, though zero concentration can also be used for particular problems.In practice, if a tracer exists only in the aqueous phase, then one can solve Eq. (4b) for aqueous diffusive transport by setting α → ∞ and using a zero-flux boundary condition at the top and bottom boundaries.

Numerical implementation
We now solve Eq. (4a) using the finite volume method.First, we discretize the equation spatially and convert it into the following ordinary differential equation ( where diag(V ) indicates a diagonal matrix formed by the vector V and Other coefficients are defined as Conductance c N+1/2 is zero when zero flux bottom boundary condition is used, but here it is included to enable Eq. ( 9) to accommodate the bottom boundary condition given as a constant tracer concentration (C b ).For this latter case, one can simply set c N+1/2 to c N −1/2 .The ODE system formed by Eq. ( 9) (together with Eq. 10) can be easily solved with various temporal discretization methods.For instance, the ODE solvers (e.g., ODE45, ODE23) in MATLAB provide a very straightforward way to obtain the solutions.In addition, we point out that by solving the equation implicitly with a very large time step (as we have done for our evaluation against steady-state analytical results below), Eq. ( 9) can also provide the steady state solution that has been used in several models to derive rates of soil methane consumption and production (Curry, 2007;Zhuang et al., 2004Zhuang et al., , 2013)).However, our formulation is more general and can be applied to model multiple gas species simultaneously.
The aqueous-primary equation Eq. ( 4b) can be solved analogously, but it should be processed with appropriate definitions of the conductances.

Evaluation with analytic models
We used two steady-state analytic models and one transient analytical model to evaluate the spatial discretization in Eqs. ( 9) and ( 10).
The first analytical model is the steady-state CO 2 diffusion model presented in Tans (1998), but we added aqueous diffusion, which was not included in his Eq. ( 2).The governing equation of the modified Tans model is whose solution is We list the parameter values used in our example application in Table 2.
We craft the second model to mimic the methane cycle in a peatland, with methane consumption in the unsaturated topsoil (defined as from the soil surface to depth z 1 , below which the soil is saturated) and a constant methane production from depth z 1 to z 2 .We could have replicated published methane models and compared with porewater CH 4 concentrations, but other uncertainties (e.g., uncertainties in parameterization, formulation, and measurement) would obfuscate a direct evaluation of the accuracy of our diffusive transport numerical formulation.
The governing equation of the steady-state methane model is whose solution is found (see Supplement) as where Parameter values used in our example application are listed in Table 3.The transient model considers the release of a tracer from a constant point source (C 0 ) into a media, which is connected to another media at some distance z 1 .The tracer has different diffusivities and solubilities in the two media, and its concentration is kept zero at the bottom of the second media.The model solves for the temporal evolution of the tracer in both media.Mathematically, the model is formulated as where we note the variables in Eq. ( 15) are now not necessarily related to gaseous and aqueous phases, but we simply keep the denotations for simplicity.The initial condition to Eq. ( 15) is set as zero tracer concentration everywhere inside the column.The model represented by Eq. ( 15) can represent a few different problems, such as contaminant diffusion from human skin into blood (e.g., Riley et al., 2004) or heat conduction between two metals of an alloy (e.g., Carslaw and Jaeger, 1986).When all coefficients are given as constant, Eq. ( 15) has the analytic solution (Carslaw and Jaeger, 1986): where s = z−z 1 , k = D g /D w .The eigenvalues, β n , are solutions of cos (βz 1 ) sin (kβL) which is solved by Newton iteration methods.In evaluating Eq. ( 16), we only used the first 17 eigenvalues, because more eigenvalues did not significantly improve the estimation.Parameter values for the example application are listed in Table 4.In all numerical experiments, we discretized the vertical soil profile using a scheme modified from CLM4.5 (Oleson et al., 2013), which defines the node depth of layer j as and the thickness of each layer as and the depth at interfaces as For all numerical experiments, the total soil column depth is set to 3.7 m.For the numerical approximation to the transient problem, the Crank-Nicholson method (Crank and Nicholson, 1947) is used for temporal discretization.

Results and discussion
Driven by the prescribed CO 2 source (Fig. 2a), the numerical solution using 100 numerical layers predicts a soil CO 2 profile very close to the exact solution throughout the soil column (Fig. 2b).The maximum relative error is about 0.2 %, which is hard to discern visually.Decreasing the number of numerical layers to 20 leads to visually discernible deviations  from the analytic solution and the maximum relative error increases to ∼ 4 %.However, the maximum relative error in both cases is near the surface, where the CO 2 concentration is low.The 20-and 100-layer simulations predict a surface CO 2 efflux of about 1 and 0.03 % accuracy, respectively, with respect to the analytical flux.These results indicate our numerical technique is sufficient for most soil dual-phase diffusion modeling applications.
The evaluation of the CH 4 numerical solution against the analytical CH 4 model in general shows good accuracy.As for CO 2 , more numerical layers lead to better numerical accuracy.Both 20-layer and 100-layer solutions indicate significant deviations from the analytical soil CH 4 profile, but the largest relative error is about 5 % for the 20-layer simulation and less than 1 % for the 100-layer simulation.The largest relative error occurs at the bottom of the topsoil (10 cm) in both cases.Because of the numerical approximation, both numerical solutions do not have numerical layers interfaced at 10 cm, which, when combined with the abrupt transition from the first-order consumption in topsoil to the constant methane production rate between 10 cm and 40 cm, lead to the largest relative error.Nevertheless, considering that the 20-layer and 100-layer solutions have 7 and 3 % relative error, respectively, in approximating the surface methane fluxes, the numerical algorithm should again satisfy the needs of modeling methane dynamics in ecosystem biogeochemical models.
For the transient model, we first compared the temporal evolution of tracer concentrations at three depths (7, 15, and 200 cm) (Fig. 4a, b, and c, respectively).Again, both the 20layer and 100-layer simulations show visually very accurate results (with mean relative error less than 5 % for the 20layer and less than 1 % for the 100-layer), though there are large errors (> 100 %) in the first 10 minutes of the simulations (when tracer concentrations are very low), which are probably caused by errors in both the finite volume approximation and the numerical error in evaluating the analytic results (see Supplement Fig. S1 for the latter case).When the steady-state solutions are compared (Fig. 4d), the two numerical models show mean relative error within 2 %, indicating our numerical algorithms are very accurate.In both transient and steady-state cases, the numerically predicted top surface fluxes agree with the analytical solution with a mean relative error within 10 % for the 20-layer and 2 % for the 100-layer simulations, respectively, after the first 20 minutes of simulation.To summarize from the three evaluations, we contend Eq. (4a) and its approximation Eq. ( 9) should be very helpful for solving the dual-phase diffusion problem.

Conclusions
Dual-phase diffusion is an important process that needs to be represented in depth-resolved biogeochemical models.Here we reviewed existing formulations in the literature and categorized three forms used in biogeochemical models.We recommend that the gas-primary form (Eq. 4a) is the most convenient to solve.Its finite volume approximation can represent tracer transport in variably saturated soils without the need for special treatment of the air-water interface (as is often done in existing methane models).Our evaluation of the numerical algorithm with three analytical models demonstrated good accuracy of our numerical model, with some dependence on spatial discretization.We hope our results can help researchers develop simple but mechanistic models for some scientific questions where more complex reactivetransport models are not necessary.
The Supplement related to this article is available online at doi:10.5194/bg-11-3721-2014-supplement.

Figure 2 .
Figure 2. Comparison between the numerical and analytical solutions for the steady-state soil CO 2 model: (a) net CO 2 source profile, as specified in the second term of Eq. (11); (b) analytical and predicted soil CO 2 profiles.N20 indicates solution using 20 numerical layers and N100 indicates solution using 100 numerical layers.At the top boundary, the CO 2 concentration equals 350 ppmv.

Figure 3 .
Figure 3.Comparison between the numerical and analytical solutions for the steady-state soil methane model: (a) net CH 4 source profiles and (b) soil gaseous CH 4 profiles.N20 indicates the solution using 20 numerical layers and N100 indicates the solution using 100 numerical layers.

Figure 4 .
Figure 4. Comparison between the numerical and analytical solutions for the point source tracer diffusion model: (a) temporal evolution of tracer concentration at 7 cm; (b) temporal evolution of tracer concentration at 15 cm; (c) temporal evolution of tracer concentration at 200 cm; and (d) steady-state solution.N20 indicates the solution using 20 numerical layers and N100 indicates the solution using 100 numerical layers.

Table 1 .
An incomplete literature survey of different formulations that have been used for modeling dual-phase diffusive transport.Bulk soil CH 4 concentration is the primary variable.Saturated and unsaturated soil use different but constant diffusivities D i .

Table 2 .
Parameters used for the steady-state CO 2 model.

Table 3 .
Parameters used for the steady-state CH 4 model.

Table 4 .
Parameters used for the transient tracer diffusion model.