Improvement of stomatal resistance and photosynthesis mechanism of Noah-MP-WDDM (v1.42) in simulation of NO2 dry deposition velocity in forests

The rapid urbanization and economic development of China has led to a dramatic increase in nitrogen oxide (NO2) emissions, causing serious atmospheric nitrogen pollution and relatively high levels of nitrogen deposition. However, despite the importance of nitrogen deposition, dry deposition processes in forested areas are still insufficiently represented in current global and regional atmospheric chemistry models, which constrains our understanding and prediction of spatial and temporal 5 patterns of nitrogen transport in forest ecosystems in South China. The offline 1-D community Noah land surface model with multi-parameterization options (Noah-MP) is coupled with the WRF-Chem dry deposition module (WDDM) and is applied to further understand and identify the key processes that affect forest canopy dry deposition. The canopy stomatal resistance mechanism and the nitrogen-limitings scheme for photosynthesis in Noah-MP-WDDM are modified to improve the simulation of reactive nitrogen oxide dry deposition velocity. This study finds that the combined improved stomatal resistance mechanism 10 and nitrogen-limitings scheme for photosynthesis (BN-23) agrees better with the observed NO2 dry deposition velocity, with mean bias reduced by 50.1%, respectively. At the same time, by comparing the different mechanisms of the two processes of canopy stoma resistance and leaf nitrogen-limiting factors, this study also finds that the diurnal changes in dry deposition velocity simulated by each regional model present four sets of distributions. This is mainly due to the different ways that each integrated mechanism handles the opening and closing of stomata at noon and the way the nitrogen-limiting factor acts. 15


Introduction
Transport and deposition of nitrogen-containing compounds is one of the most critical processes in the study of biogeochemical cycles (Gruber and Galloway, 2008). Atmospheric nitrogen deposition is not only the main way that atmospheric reactive nitrogen is removed, but is also an important source of nitrogen for ecosystems (Jefferies and Maron, 1997;Horii et al., 2005).
Nitrogen deposition affects changes in the carbon sink in forest ecosystems by affecting plant growth and death (De Vries et al., 2009;Bernhard, 2012). An increase in nitrogen deposition will cause an increase in litter and a decrease in soil decomposition, which will increase the carbon fixation of the soil (Stevens et al., 2004;Liang et al., 2020). Meanwhile, soil acidification caused by nitrogen deposition will reduce the number of microorganisms in the soil, reduce the production of methane, cause the degradation of peatland, and jointly affect the balance of greenhouse gases and the climate (Xu et al., 2009;Cui et al., 2010; 25 Seinfeld and Pandis, 2012;Erisman et al., 2014). At present, studies have shown that there has been a sharp rise in global and regional atmospheric nitrogen deposition, which exceeds the critical load of local ecosystems in many regions (Liu et al., 2013;Yu et al., 2019).
In order to evaluate the impact of atmospheric nitrogen dry deposition on ecosystems, it is important to accurately estimate the dry deposition fluxes of nitrogen components (Wu et al., 2011(Wu et al., , 2012Tian et al., 2018). Scholars calculate the dry deposi-30 tion velocity of nitrogen-containing components or estimate the effect of variation on dry deposition flux based on the global or regional numerical models (Phillips et al., 2006;Zhao et al., 2017;Han et al., 2017;Zhong et al., 2020). The biased results of nitrogen deposition from modelling compared to those from observations range from -70% to 800% (Chang et al., 2020a).
Part of the estimated deviation comes from the nitrogen emission inventory input bias in the model simulation. For example, Galloway et al. (1994) predicted the nitrogen deposition pattern for 2020, and a large deflection area appeared in an area where 35 emissions did not increase as expected. Another part of the deviation comes from the inaccurate simulation of nitrogen concentration. In this situation, the simulated concentration can be verified and nudged by measuring stable oxygen and nitrogen isotope ratios (Guerrieri et al., 2020). At the same time, the simplification and biases from the deposition mechanism in the models can not be ignored compared with satellite retrievals . Deposition velocity is difficult to measure and is affected by many coupled physical, chemical and biological processes occurring at the deposition interface. 40 Therefore, the resistance-velocity method, which is similar to Ohm's law, is used to calculate the dry deposition velocity of various species between the surface and the atmosphere (Szinyei, 2015). In this method, the dry deposition velocity (V d ) of gaseous matter is expressed as the reciprocal of the total resistance (R t ) of the process of atmospheric pollutant deposition to the surface. The total resistance is determined by aerodynamic resistance(R a ), quasi-laminar boundary layer resistance (R b ), and canopy resistance (R c ). Their relationship is generally characterized by Equ. 1.
Here, R a is calculated by the micrometeorological parameters, which depend mainly on local atmospheric turbulence intensity while R b is driven by the diffusion coefficient and air viscosity of gaseous matter. The calculation of these two resistances is basically similar in different deposition resistance mechanisms (Finnigan, 2000). At present, the treatment of dry deposition processes affected by turbulent diffusion in numerical models includes two parts: one is the turbulent diffusion process from 50 the bottom of the atmospheric boundary layer to the canopy, and the other is the turbulent exchange process inside the canopy (Flechard et al., 2011(Flechard et al., , 2013. The calculation of R a is usually based on the turbulent transport part of the land surface model. Most current models are based on the near-surface-layer similarity theory, which first calculates the surface roughness and zero plane displacement, and then calculates the turbulent transport coefficient according to the flux gradient relationship under different stratifications 55 (Makar et al., 2017). The calculation of turbulent exchange inside the canopy is more complex, and is highly related to the structure of the vegetation canopy and other local properties (Finnigan et al., 2009). Some forest fire models are based on the measured empirical wind speed profiles in the canopy, and other models use the assumption of neutral stratification to solve the turbulent flow fields of the canopy, such as SSiB, SVAT, and BATS, etc (Yongjiu and Qingcun, 1997;Yang and Friedl, 2003;Falge et al., 2005;Moon et al., 2019).

60
Furthermore, the calculation of R c is more complicated and diverse than that of R a , because R c is closely related to difference in the underlying surface, vegetation, soil, and other conditions (Wu et al., 2018). Due to different considerations of the underlying surface, R c is usually further decomposed based on canopy structure, surface properties of deposition receptors, biochemical reactions of deposition materials, and other canopy processes (Massad et al., 2020). For the surface of the vegetation canopy, models are refined to consider the resistance of the stomata, mesophyll, epidermis, soil, and other canopy surface factors (Dai et al., 2004;Massad et al., 2020). For example, the multi-layer forest canopy model is used to calculate the canopy stomatal resistance layer by layer at monitoring sites in the Clean Air Status and Trends Network (CASTNET) (Li et al., 2016). As a counterexample, the ocean, which was thought to be a relatively simple surface, has evolved from consideration of smooth levels to sea surface fragmentation, different particle humidities, and other factors (Schulz et al., 2012). However, our current understanding of the exchange of nitrogen oxides between the atmosphere and biosphere remains incomplete (Delaria 70 and Cohen, 2020). For instance, the possible existence of a NO 2 compensation point toward the leaf surface in forests has been controversial as a result of experimental comparison . The coupling of canopy photosynthesis, nutrient stress, the impact of mesophyllic processes, and other plant physiological processes is still poorly resolved in the field of dry deposition model improvement (Massad et al., 2020).
In this study, we apply different improved stomatal resistance mechanisms and nitrogen limitations on photosynthesis mech-75 anisms to the Noah-MP model coupled with dry deposition schemes to explore changes in nutrient stress in stomatal conductance and evaluate the consequences of these changes on NO 2 dry deposition velocity. This paper is organized as follows: besides this introductory Sect. 1, Sect. 2 presents a full description of the improved stomatal resistance mechanisms and different schemes of nitrogen limitation of photosynthesis. Sect. 3 includes model evaluation and discussion about the influence on NO 2 dry deposition velocity simulation, the respective path of canopy stomatal and photosynthesis processes and sensitivity 1, 3, 3, 2, and 2 in the above physical parameterization schemes (Chang et al., 2020b).

Coupling of Stomatal Resistance Schemes
Previous studies have generally used the Jarvis stomatal conductance model, which is based on environmental factors such as photosynthetic effective radiation, temperature, humidity, and soil water to calculate canopy stomatal resistance (Jarvis,100 1976). Compared with Jarvis, the Ball-Berry stomatal conductance model (Ball et al., 1987) is a better mechanism for the stomatal resistance of the canopy; it calculates the stomatal resistance based on the through-canopy photosynthesis rate, CO 2 concentration, and humidity on the leaf surface as shown in (2). This type of mechanism requires a coupled photosynthesis model to calculate or observe the photosynthesis rate of the canopy, and the photosynthesis model depends on the setting of many key plant physiological parameters (such as optimal photosynthesis efficiency, catalytic enzyme activity parameter 105 Q10, etc.). The Noah-MP model sets the stomatal conductance slope of the Ball-Berry mechanism as a constant, which is not suitable and will cause a large simulation bias. Therefore, we integrate observational experimental results, statistical fitting or plant physiological model equations in photosynthesis, stomatal conductance and other aspects of plant physiology in this study, writing the equation as a subroutine in the coupled single-point model.
The calculation equation is Ball et al. (1987) as follows: where R ac is the in-canopy aerodynamic resistance, which is common to all gases and R g and R cut are the resistances for the uptake by the ground/soil and canopy cuticle. Similar to the work of Wesely (1989), R g and R cut are parameterized for O 3 from look-up tables.
The equations integrated into the single-point mechanism model are shown in Table 1 - Leuning (1990) introduced a CO 2 compensation point Γ to improve the Ball-Berry equation so that it can simulate the net photosynthetic rate and stomatal conductance when the CO 2 concentration on the blade surface is equal to the 125 compensation point (Table 1, MBM-2). Later, the method of Lohammar et al. (1980) was adopted to replace RH with the water vapor saturation function f (D) ( Table 1, MBM-4). These equations have been applied to a variety of plant physiological models (Leuning, 1995); -Aphalo and Jarvis (1993) separated the effects of temperature and water vapor difference D, which more directly reflected the effect of temperature on stomatal conductance than the original Ball-Berry equation (Table 1,  The combination of the DVEG mechanism and the Ball-Berry model can comprehensively consider the interaction between photosynthesis rate and canopy stomatal impedance. This is physiologically significant in that it balances the supply and demand of CO 2 in the chemical reaction of photosynthesis, so as to maintain a reasonable concentration of CO 2 in the mesophyll tissue. However, at the vegetation canopy scale, the photosynthetic rate is also related to the nitrogen content of leaves. Currently, the commonly used biogeochemical models usually express the effect of nitrogen on the photosynthetic 145 process based on the relevant theory of nitrogen limitation, but they often simplify it (Li et al., 2013).
In this study, the original DVEG mechanism of Noah-MP set the nitrogen limitation factor of leaves f (N ) as a function of leaf nitrogen concentration (C leaf N ) and the maximum nitrogen concentration parameters (FOLNMX) of this vegetation type (f (N ) = C leaf N · F OLN M X −1 ). However, C leaf N and FOLNMX were set as two constants which is obviously an oversimplification for land surface simulation in a large area (Bonan, 1995). For different types of plants, the nitrogen content in 150 leaves should make a significant difference in photosynthetic nitrogen utilization efficiency (Zheng and Shangguan, 2007). For regional nitrogen deposition simulation, it is obviously inappropriate to simplify the nitrogen-limiting process in leaves, so a more accurate description of the effect of nitrogen on plant photosynthesis and a more accurate estimation of the effect of nitrogen deposition on the whole forest ecosystem are needed.
According to whether the nitrogen content of plant tissues is directly taken as the variable in the equation, the current 155 expressions of how nitrogen affects photosynthesis (as shown in Table 2) can be divided into implicit and explicit expressions: -Implicit For the photosynthetic rate model calculated by the Farquhar Model (Farquhar et al., 1980), the photosynthetic rate is determined by the minimum value of carboxylation efficiency (W c ), carboxylation efficiency (W j ) and organophosphorus carboxylation efficiency (W e ), which is limited by the concentration of chlorophyll photoenzyme (Rubisco), in which 160 W c and W e are proportional to the maximum carboxylation rate (V cmax ).
Therefore, the effect of nitrogen on photosynthesis, V cmax , is reflected mainly in the limitation of f (N ). As mentioned above, f (N ) was set by two constants in the DVEG dynamic vegetation process mechanism. In addition, models such as AVIM (Ji, 1995), CLM4.0 (Oleson et al., 2010) and Noah-LSM (Bonan, 1995) also directly take f (N ) as a parameter, ranging from 0.5 to 1( Table 2,   However, BEPS (Liu et al., 1999), DLEM (Tian et al., 2010), IBIS (Liu et al., 2005), InTEC (Chen et al., 2000) and other models use the ratio of the optimal carbon-nitrogen ratio (B V max) to the simulated actual carbon-nitrogen ratio (B L ) to represent f (N ) ( Table 2, MNM-2).

185
After integrating all the improved equations of the canopy stomatal resistance mechanism and the nitrogen-limiting schemes for photosynthesis into the single-point model in the form of subroutines, an orthogonal experimental scheme (Table 3) was adopted to simulate them, and all the experimental schemes were driven by the same meteorological forcing data. Since all the mechanisms can be combined into 42 combinations, the current version number is set at v1.42 in this study.  Table 4. It can be seen that the simulation of average SH is overestimated by about 20 W · m 2 while the average LH is underestimated by about 0.1 W · m 2 compared with observations. The models perform reasonably 195 well for most simulations with uncertainties within a factor of 0.5 ∼ 2 (Fig. 1).

Performance of Vd Simulation with Different Mechanisms
The model simulation show obvious underestimation of V d . The simulated average V d is about a quarter of the observed results.
And the correlation coefficient is very low and basically cannot reflect the trend characteristics (Fig. 2). On the one hand, the Noah-MP-WDDM model itself has a poor ability to simulate the change trend of deposition. On the other hand, it is also 200 affected by too much precipitation in subtropical regions, poor quality control of dry deposition observation data and many missing values .
It can be seen from Fig. 2 that the simulation effect of each model mechanism is relatively poor, especially for all the combinations corresponding to the MBM-5, MBM-6, MBM-7, and MNM-5 series; the simulated V d in these series is basically concentrated around 0.05 cms. This indicates that the stability of the parameterization of these series of mechanisms is relatively 205 high, and the disturbance caused by different schemes in other processes is suppressed.
There is a magnitude difference between the results of the simulation and the observed V d , which may be because these mechanisms are not supported with some coniferous species because conifers have little direct stomatal response to elevated CO 2 (Medlyn et al., 2011;Katul et al., 2012). However, from the statistical results (Fig. 3), we can still see the variation of the simulation bias caused by different mechanisms. Most of the simulated combinations underestimate the average dry deposition velocity, but only three mechanism combinations-BN-16, BN-26, and BN-36-overestimate it. It can be seen that the average simulation deviation of BN-23 is the lowest among all mechanism combinations. Compared with the default BN-11 mechanism, its average deviation is reduced 215 from -0.0371cm/s to -0.0185cm/s, which reduces the relative deviation about 50.1%.

Implications for NO 2 Dry Deposition Velocity Diurnal Simulation
Although the trend simulation ability of the models is statistically weak, and the absolute difference in the average dry deposition velocity obtained from simultaneous results is small, we can still see that the model captures certain dry deposition characteristics from the daily cycle changes. Fig. 4 uses a daily variation curve to show the simulation results of the effects of 220 each mechanism combination on the dry deposition velocity. It can be seen that for the daily variation of N O 2 dry deposition velocity, different mechanisms still show considerable pattern differences. Fig. 4 is the daily change in the observed value. Note the large fluctuation of its standard deviation, indicating a large fluctuation during the V d observation period. This is because of the turbulent exchange caused by the fragmentation of the boundary layer inside and outside the mountain forest canopy and the effect that the change in atmospheric stability has on 225 the turbulence . The black line in Fig. 4 is BN-11 and the green line is BN-23. It can be seen that compared to the original BN-11, BN-23 shifts the V d simulation upward mainly during the day. At the same time, it can be seen that the standard deviation of the observations basically covers all the standard deviations of the simulation results. This reflects the stability of the models, which means that these improved mechanisms should show similar performance when transplanted to other, similar types of underlying surfaces.

230
In addition, it is worth noting that some deposition observation studies believe that the V d value at midday is the most noteworthy (Kavassalis and Murphy, 2017;Ke et al., 2020). Therefore, we also pay attention to the mechanism with the smallest simulation deviation at midday, which is BN-46. It can be seen that the simulated value of BN-46 is basically consistent with the V d observation, with a bias of about 0.001cm/s. (represented by BN-23), and the accurate-at-noon group (represented by BN-46). The original Noah-MP-WDDM  belongs to the same group as BN-23 because their theoretical assumptions are consistent. The appearance of this grouping is quite interesting, because it illustrates that there are relative differences in theoretical assumptions about stomatal resistance and nitrogen limits for photosynthesis. Therefore, in the next section we will try from the perspective of canopy deposition 240 resistances with the four representative combinations and the Noah-MP-WDDM default combination (BN-11), to discuss the effects of these improved scheme groups.

Aerodynamic and Quasi-laminar Boundary Layer Resistance
Since the improved mechanisms are concentrated in the canopy layer, disturbances to aerodynamic resistance (R a ) and quasi- In addition, it is worth noting that the diurnal variation of R b is more consistent with the observed V d , which echoes the hypothesis of turbulence exchange caused by the breakage of the inner and outer boundary layers of the mountain forest canopy and changes in atmospheric stability proposed by Zhang et al. (2017). For this reason, we hope that in the future, the model will be able to express these situations by accurately expressing the forest structure.

Canopy Resistance
The difference in the simulation of the canopy resistance (R c ) of each improvement scheme is the main source of the difference in the simulations of the dry deposition process, as shown in Fig. 7. It can be seen that for the consistently underestimated group represented by BN-55, the underestimation of deposition velocity comes from a large overestimation of R c , indicating that the assumptions of these mechanisms are not suitable for subtropical forests as a whole. It is possible to draw inappropriate 260 conclusions on such underlying surfaces from the source model results.
For it can be seen that their theoretical hypotheses are the same, which can effectively reflect the physiological process of the increase in stomatal resistance caused by the closure of stomata at noon. The degree of stomata reopening in the afternoon is slightly weaker, which makes the diurnal dry deposition velocity curve high in the morning and low in the afternoon. This is correct in a general sense, but there is a certain mismatch between the simulation and the observed results

265
(low in the morning and high in the afternoon). We estimate that because the sample site where the flux tower of Dinghushan is located on a westward slope, the physiological activity of the vegetation canopy is weaker in the morning and stronger in the afternoon. This indicates that for model improvement, the parameterization of the difference between sunlit and shaded leaves should be strengthened; otherwise it will be difficult to express this phenomenon.
It can be seen that the theoretical hypotheses for BN-26 and BN-46 are also the same, but neither can reflect the closure of 270 stomata at noon. The difference is reflected mainly in the intensity of the decrease in R c during the day, and the amplitude of the disturbance to the deposition velocity is greatly enhanced when the deposition resistance is lower than 1000 s/m. It can be seen from Fig. 4 and Fig. 7 that the difference in the deposition velocity curves obtained by the simulation of  is not yet apparent at 07:00 in the morning. At 12:00, the difference in canopy resistance, which is only about 200 s/m lower, causes the BN-26 group to greatly overestimate deposition velocity. Under this series of mechanisms, the misestimation 275 or disturbance of key parameters is likely to change the expected results.

Conclusions
In using Noah-MP-WDDM to study dry deposition processes, we implemented new features and applied several corrections to the code. Compared to Noah-MP-WDDM v1, the improvement of the canopy stomatal resistance mechanism and the nitrogenlimiting schemes in Noah-MP-WDDM v1.42 result in much improved agreement with measurements from nitrogen dry depo-280 sition results. Our tests for Dinghushan show that Noah-MP-WDDM v1.42 now gives a better deposition velocity simulation due to modification of the canopy stomatal resistance mechanism and the nitrogen-limiting schemes for photosynthesis.
Our results emphasize the importance of the canopy stomatal carbon dioxide compensation mechanism (Leuning, 1990) and the GPP-controlled leaf nitrogen-limiting factor (Lin et al., 2000) for the simulation of nitrogen deposition. Considering the combination of these two mechanisms (BN-23 schemes in Noah-MP-WDDM v1.42 instead of Noah-MP-WDDM v1), the 285 average simulation bias is reduced about 50.1%.
Our discussion shows that the canopy stomata and leaf nitrogen-limiting mechanisms from various classic models cannot well express the diurnal changes in leaf canopy resistance, especially the underestimation in the daytime, by the combination of the (Yu et al., 2004) and (Thornton et al., 2002) mechanisms (BN-55) and the effect of the (Cao and Woodward, 1998) mechanism on stomatal closure (MNM-6) at noon. This may be a source of bias in the simulation of nitrogen deposition flux 290 by these mechanisms' source models.