Tracking the moisture sources of an extreme precipitation event in Shandong, China in July 2007: A computational analysis

This paper utilizes a modified Water Accounting Model (WAM) to track the moisture sources of an extreme precipitation event in Shandong during 18–20 July 2007. It is found that different methods in dealing with the residual of the water budget always produce different results in moisture recycling calculations. In addition, results from the backward tracking without the residual are in complete agreement with those from the forward tracking with the residual, and vice versa, implying a mathematical consistency. We thus analyze and derive the conditions under which the two tracking approaches equate with each other. We applied the backward tracking to the Shandong extreme rainfall case and obtained quantitative estimates of moisture contributions of three selected regions away from the rainfall area. The results indicate that the spatial pattern rather than numerical value of the recycling moisture is more reliable in tracking the moisture sources. The moisture of this Shandong rainfall event comes mostly from the nearby upwind area in Southwest China, which is of the terrestrial origin; while the moisture originating from the neighboring West Pacific contributes little to this event.


Introduction
Research into the moisture sources of precipitation is important for understanding water cycling (Ralph et al., 2006). Such research reveals the mechanisms by which water is evaporated from one area, transported, and falls in another area as precipitation. This process is termed moisture recycling or precipitation recycling. There are four categories of moisture recycling models: bulk models, general circulation models (GCMs) with tagged water, Lagrangian methods, and isotopic analysis.
GCM with tagged water is a tagging technique in the Eulerian frame. GCMs implement the numerical water vapor tracers (WVTs) that experience the same processes as atmospheric water (Bosilovich and Schubert, 2002;Gimeno et al., 2012). The evaporation, transportation, and sink destination of WVTs are tracked and recorded. GCMs may be the most comprehensive and complex models in atmospheric science; they provide the most physical meaning and require the most complex computation. The accuracy of GCMs depends on the simulation or parameterization of various atmospheric processes, such as convection, turbulence, and longwave and shortwave radiation.
The Lagrangian-type methods are divided into two kinds. One is the back-trajectory method, which considers water vapor as a passive tracer along quasiisentropic surfaces (Dirmeyer and Brubaker, 1999;Dirmeyer and Brubaker, 2007). This method does have some limitations. Kurita et al. (2004) noted that the adiabatic assumption of this method may not be valid in summer over land. Moreover, precipitation is divided into many parcels of the same mass that originate from an atmospheric level according to a probability distribution model, which may not be realistic (Fitzmaurice, 2007). The other kind of Lagrangian Supported by the National Science and Technology Support Program of China (2012BAC20B06). * Corresponding author: zhangchi.abc@gmail.com. ©The Chinese Meteorological Society and Springer-Verlag Berlin Heidelberg 2014 methods is based on atmospheric particle dispersion models, which have been incorporated into tools such as FLEXPART (FLEXible PARTicle dispersion model; Stohl et al., 2005) and HYSPLIT (HYbrid Single Particle Lagrangian Integrated Trajectory Model; Draxler and Hess, 1998). The atmospheric particle dispersion models also have limitations. They can only account for the net value of "evaporation minus precipitation" in air parcels, which means that they cannot separate evaporation and precipitation. This causes bias in the source-receptor relationship (Gimeno et al., 2012;van der Ent et al., 2013).
Isotopic analysis is a method that can validate recycling models with real data. The studies by Kurita et al. (2003Kurita et al. ( , 2004 and associated circulation model studies by Yoshimura et al. (2003Yoshimura et al. ( , 2008 provided quantitative information on evaporative sources of precipitation. However, the isotopic data are not yet sufficient for this method, and there are large modeling uncertainties such as fractionation parameters associated with the isotopic analysis (Gimeno et al., 2012).
Bulk models are well known for their computational simplicity and flexibility in defining regions (Bosilovich and Chern, 2006). They are divided into two categories: analytical models and numerical models. Budyko and Drozdov (1953) proposed the first analytical model, i.e., the Budyko model, which is the basis of various bulk models. Since then, many bulk models have sprung up, such as those of Grigor'eva (1965), Brubaker et al. (1993), and Burde and Zangvil (2001). Each model deals with different situations and has its own limitations. The early bulk models neglect the change in atmospheric moisture storage and consider it to be small compared with other terms on longer timescales, thus they are not suitable for moisture recycling on shorter timescales. Dominguez et al. (2006) inserted the moisture storage term back into the basic water conservation equation and developed the dynamic recycling model (DRM), in which the conservation equation is analytically solved under a Lagrangian frame. DRM can be used on short timescales that are longer than the boundary layer mixing time (Dominguez et al., 2006). Eltahir and Bras (1994) proposed the first grid-based numerical method on moisture recycling. The basic assumption is that the ratio of local to total moisture out-flux in a grid box is equal to the ratio of local to total precipitation in the same grid box. Numerical models address the moisture recycling results of spatial variation. Later, Kurita et al. (2003), Fu et al. (2006, and Fitzmaurice (2007) incorporated the storage term into the Eltahir and Bras model, allowing it to deal with moisture recycling on the sub-monthly timescale. These models only provide the local water sources; they do not include the remote water sources for a region (Bosilovich and Schubert, 2002). Modern numerical methods overcome this. They are commonly based on the basic atmospheric moisture balance equation and the well-mixed assumption, which states that the advected and local evaporated moistures are well-mixed so that each water molecule has the same probability to be precipitated out. This is a basic assumption of bulk models and also a most questionable one that always results in a lot of discussion (Burde, 2006;Goessling and Reick, 2013;van der Ent et al., 2013). The colored moisture analysis (CMA) of Yoshimura et al. (2004), the water accounting model of van der Ent et al. (2010) and van der Ent and Savenije (2011) are representatives of modern models. These modern models can conduct not only the tracking forward but also the tracking backward in time by making precipitation the source term and evaporation the sink term (van der Ent and Savenije, 2013;van der Ent et al., 2013). Therefore, they can be used to track the moisture sources of a precipitation event. Furthermore, the contribution of the source to the sink can also be calculated by forward tracking. In theory, the contribution from source to sink calculated by the forward tracking and that by the backward tracking should be equal. However, in practice, when the data are not accurate and the water balance equation does not close, will the results still match? What are possible biases in numerical computations? In this paper, we try to answer these questions through numerical experiments and theoretical derivation. An extreme precipitation event that occurred from 18 to 20 July 2007 in Shandong Province is used as a case study. A modified water accounting model (WAM) of van der Ent et al. (2010) and van der Ent and Savenije (2011) is adopted, and the moisture sources for this Shandong precipitation event are investigated.
This paper is organized as follows. Section 2 describes the data and the precipitation event. Section 3 provides the mathematical frame and numerical implementation for forward and backward moisture tracking methods. Section 4 presents the results. Section 5 is an analysis and discussion of the methods, errors, and results. Section 6 provides the conclusions.

Data
In this study, most meteorological data were taken from the ERA-Interim (ERA-I) reanalysis at a grid of 1.5 • latitude × 1.5 • longitude from 1 to 20 July 2007. The data were downloaded from http://dataportal.ecmwf.int/data/d/interim − full − daily/levtype =pl/. These data include specific humidity, zonal and meridional wind speeds at the 23 lowest pressure levels (200-1000 hPa), and surface pressure. All are instantaneous values obtained at 6-h intervals. The ERA-I precipitation and evaporation are also included, which are accumulated values at 3-h intervals. A Chinese province-level division map of shapefile format is provided to delineate the study area. According to meteorological station records (Chen et al., 2011), there was widespread heavy precipitation over the entire Shandong Province from 18 to 20 July 2007. The accumulated precipitation over China during 18-20 July 2007 from ERA-I is shown in Fig. 1. As can be seen, heavy precipitation occurred in the middle and southwest of China. Precipitation in Shandong did not appear to be very heavy, which suggests that the ERA-I precipitation may not be accurate. The Global Precipitation Climatology Project 1-Degree Daily (GPCP 1DD) product for the corresponding period was used as a supplementary precipitation dataset. Since the GPCP 1DD product is of a 1degree resolution, the dataset was re-sampled into a 1.5-degree grid. The GPCP data on accumulated precipitation for 18-20 July 2007 better present the rain event in Shandong (figure omitted). Shandong Province is approximated by eight grid cells, as indicated in Fig

Equations for forward moisture tracking
Bulk models take the entirety of the water column. For an integral water column, the atmospheric water balance equation is given as follows: where where w is the precipitable water contained in the unit area column of air; E is the evaporation; P is the precipitation; q is the specific humidity; u and v stand for zonal and meridional wind velocities, and u p and v p stand for wind velocities at different pressure levels; g is the gravitational acceleration; and p s is the surface pressure. When tagging water from a specific region, according to the well-mixed assumption, there is where the subscript m refers to the tagged water from the source region. The balance equation of the atmospheric moisture from the source region is expressed as follows: When in the source region, δ = 1; otherwise, δ = 0.

Equations for backward moisture tracking
Section 3.1 provides a mathematical solution of water from evaporation to precipitation. In the sourcereceptor process, the tagged water conserves to Eq. (6). Source A's contribution to receptor B's precipitation in a period can be measured. If we inverse the time axis, when precipitation returns into the atmosphere as water vapor and evaporation sinks to the ground as ground water, the precipitated water will return to its original source(s). A similar process applies in the forward method. If we change precipitation into source and evaporation into sink term, the new atmospheric water balance equation is as follows: where u and v will also change sign. The well-mixed assumption still holds. Therefore, there is where the subscript m refers to the tagged water from the precipitation region. When precipitation reverses back to evaporative source(s), part of it sinks to the ground in proportion to its moisture fraction, which is E × w m /w. According to the basic conservation equation, i.e., Eq. (7), the atmospheric water balance equation of precipitation reversing is given as follows: where δ = 1 for the tagged source region (precipitation as source), and δ = 0 otherwise.

Numerical implementation
The model used here is a modification of WAM. The numerical implementation is similar. The numerical scheme for the moisture balance equations is an explicit one, and the advection is done with central difference scheme at the grid boundaries. The difference scheme is not a general forward-time centralspace (FTCS) method, since the precipitable water serves as a known variable, and we always substitute the calculated value at the next time step. The grid is 1.5 • × 1.5 • and approximated by trapezoids. The horizontal moisture fluxes, wu and wv, are transformed into unit cubic meter of liquid water per time step over grid boundaries; vertical fluxes, E and P , are transformed into unit cubic meter per time step over a grid cell; and precipitable water, w, is transformed into unit cubic meter over a grid cell. To calculate the column moisture and horizontal moisture flux, the specific humidity and wind speed at ground level are interpolated or extrapolated linearly.
To be in harmony with the explicit scheme, the scheme needs to fulfill the Courant-Friedrichs-Lewy condition to keep the computation stable. In this case, the condition is: |wu| + |wv| < w. The time step is set to 0.5 h in van der Ent et al. (2010). When the same time step is set in this study, the result still diverges. Diverged grids always occur in polar regions, where w is much smaller than in midlatitude grids. Therefore, in the present study, we made a modification by combining 10 rows of grids on each polar edge into one row of grids. In WAM, the north and south edges are generally left alone, since including them has barely any effect. When the north and south boundary grids are combined, these new boundaries cannot be left alone. They are put in the recycling computation in the modified WAM. Region A's contribution to region B's precipitation in the forward tracking method is calculated as follows: where P (t, x, y|B) is the precipitation over region B, w m (t, x, y|B) is the air moisture of origin A, w m (t, x, y|B) is the total air moisture, and t spans the precipitation period.
In the backward tracking method, the spatial contribution map for region B's precipitation is easily obtained. For each grid, its contribution quantity is where t spans the whole pre-precipitation evaporation time, which may begin more than 10 days before precipitation. For region A, it is The residence time for atmospheric moisture varies with time and space (Trenberth, 1998). If we suppose that the residence time for A is k days at a certain time, any moisture from A will have no influence or will be so small that is has a negligible influence on precipitation k days later. Then, if t in Eqs. (10) and (11) approaches infinity, the results of P A2B and P AconB should approach each other. In practice, the moisture always decays in an e-folding mode, and t is finite, so the two results may be approximate but not identical. Even when there are errors in the data, the basic moisture balance equation does not close, and the methods of dealing with residual will also influence the results. The residual α is added to the moisture balance equation as in Dominguez et al. (2006): and Eq. (6) changes into and Eq. (9) changes similarily. There are several ways to deal with the residual, such as considering it as a part of precipitation (Goessling and Reick, 2011) or assuming that α m /α = w m /w as in Yoshimura et al. (2004) or van der Ent et al. (2010). In this paper, the latter approach was adopted. Forward and backward difference schemes based on the simplest one-variable cases usually do not match. To reduce the imbalance associated with different difference schemes applied to unclosed data, a set of experiments were performed, where we added the residual term α in the forward tracking model but not in the backward tracking model (α m = 0 in the latter model). This is another modification we made on the WAM. Accordingly, another set of backward tracking with the residual was done for comparison. Numaguti (1999) found that the average time of water vapor residing in the atmosphere was 10 days. Other studies such as Trenberth (1998) and van der Ent and Savenije (2011) investigated extensively the residence time of water vapor in the atmosphere. They revealed that the residence time changes with both geographic location and season. In this paper, we assume that the residence time is around 20 days. In the present case, the precipitation period was from 18 to 20 July 2007. Based on the 20-day residence time assumption, the moisture before 1 July had a negligible influence on the precipitation. The evaporative source was traced back to 1 July. Furthermore, in forward tracking, the beginning date was set on 1 July and the ending date on 20 July. Then, the precipitation contribution from different sources during the extreme precipitation event was calculated. The source regions were determined according to the backward tracking described in Section 4. We also used the GPCP data in comparison with the ERA-I data.

Results
According to the ERA-I precipitation data, the accumulated precipitation in Shandong during the three days of the extreme event was 1.45 × 10 10 m 3 (58.7 mm in depth). According to GPCP data, it was 1.56 × 10 10 m 3 (63.2 mm). Through the backward moisture tracking model, the contributed moisture from each grid in every time step was recorded. The accumulated contribution to Shandong's precipitation from each grid from 1 to 20 July, based on backward tracking without the residual, is shown in Fig. 2. The source grids contribute 1.42 × 10 10 m 3 of water in total, which accounts for 98% of the precipitation, with around 0.03 × 10 10 m 3 left hanging in the air "to be precipitated." Backward tracking with GPCP precipitation was also performed. The moisture contribution pattern of GPCP is only slightly different than that of the ERA-I, with the exception that the range of the same contribution level expanded. In backward tracking considering the residual, the grids contribute 1.14 × 10 10 m 3 of water in total, with 0.02 × 10 10 m 3 of water left in the air. When the residual was added, the tagged water from Shandong precipitation did not conserve. In terms of water balance, backward tracking without the residual provides a better result.
The spatial contributions based on ERA-I with residual are shown in Fig. 3. The spatial pattern differs a little from the backtracking without residual. According to the tracking results, some water sources even traversed the Arabian Sea to reach the east coast of Africa and cross the equator. But the further they are from the sink, the smaller the contribution of these grids to Shandong precipitation, since much of the moisture is depleted on the way. The higher contribution values appear upwind near the precipitation region, where the local evaporation and the tagged moisture ratio are high. The contributions are more concentrated on land, especially from northern to southwestern China.
Since the results of the forward tracking method must be somewhat consistent with those of the back-  ward tracking method, three sample squares were picked as water sources to perform forward tracking. These regions were picked according to the grid contribution gradient. They belong to three different contribution levels, as shown in Fig. 2. Each comprises 4 × 4 grid cells. According to the results of backward tracking shown in Fig. 2, the three regions contribute (from near to far) about 1.92 × 10 9 , 1.06 × 10 9 , and 3.53 × 10 8 m 3 of water respectively, accounting for 13.3%, 7.3%, and 2.4%, respectively.
For GPCP data with the same residual scheme, the resulting contributions are 1.93 × 10 9 , 1.12 × 10 9 , and 3.74 × 10 8 m 3 (12.3%, 7.2%, and 2.4%), respectively. For the ERA-I data considering the residual, the results are 1.71 × 10 9 , 8.67 × 10 8 , and 1.75 × 10 8 m 3 . As mentioned in previous sections, in order to reduce the imbalance effect of different difference schemes on imbalanced data, the forward tracking method considering the residual was implemented for comparison with the backward tracking method without the residual. Applying the forward tracking model, the three source regions' contributions are 1.92 × 10 9 , 1.06 × 10 9 , and 3.53 × 10 8 m 3 , which are the same as the results calculated by the backward track-ing method. The GPCP forward tracking experiment produce results of 1.93 × 10 9 , 1.12 × 10 9 , and 3.74 × 10 8 m 3 . The two results agree decently.
Region 2's contribution based on ERA-I precipitation data from 18 to 20 July is used as an example in Fig. 4. The evaporated water leaves the source region and precipitates wherever it flows. As seen in Fig. 4, most of region 2's evaporation does not fall into its source. When there is heavy rain in the downwind direction and the tagged moisture ratio is high, the recycled precipitation from region 2 is also high. The major sink is in the middle of China, downwind of the source region, where precipitation is the highest (ERA-I dataset; Fig. 1) and gradually reduces as moisture flows away.

Analysis and discussion
It does not seem to be coincidental that the results of forward and backward tracking match completely in the Shandong rainfall case. Several other experiments conducted by us also support the idea that forward tracking with the residual agrees with backward tracking without the residual. There seems to be a mathe- matical necessity in the difference scheme. According to experiments conducted on this issue, we could express this relationship accurately as follows: Source A's evaporation at time t 1 has an influence (contribution) f 1 (through forward tracking method with residual) on receptor B's precipitation at time t 2 . Then, through backtracking without the residual, the precipitation of B at time t 2 reverses its way back to source A. The contribution of source A at time t 1 denoted as f 2 equates with f 1 .
The analysis of such situations is complex. When moisture from source A spreads to sink B, there are many paths of A reaching B. We attempt to simplify this complexity and consider just one of the paths. Other cases are additions to or derivations of this situation. Meanwhile, the flows of the grids along the path are various. To further simplify, consider that there are only flows along this path, and the flows are in a single direction from A to B. Then, the question becomes how the evaporation E of A at time t 1 influences the precipitation P of B at time t 2 in this flow path. However, it is still hard to make direct comparisons between forward tracking and backward tracking results, since the iterations are very complex. Within several time steps, the iterations will make the expressions extraordinarily large, and it will be difficult to compare the expressions directly. We try to sort out the analysis through the following steps.
Step 1: The grid next to A is denoted as N 1 . Prove that the statement holds between grid A and its neighboring grid N 1 , which is much easier to obtain. Then the statement holds between grid N 1 and the next grid N 2 .
Step 2: Prove that the statement holds between grid A and grid N 2 . A's evaporation transfers a volume to grid N 1 , denoted as E 1 , and N 1 and N 2 conform to the statement. E 1 can be viewed as evaporation from N 1 . In the next time step, the left evaporation of E in A transfers another volume to N 1 , which is also viewed as evaporation of N 1 . Then, the statement between A and N 2 can be precisely deduced, so is N 1 to N 3 , N 2 to N 4 , . . ., and so on.
Step 3: If the statement between N 1 and grid N i stands, it is inferred that it stands between A and N i through Step 2 with N 1 being a pivot point. Then, along the path from A to B, the statement stands be-tween A and B. This completes the demonstration of the steps. In the demonstration, there is a stronger condition than stability to meet when the two results completely agree; that is, |wu| + |wv| + E < w. In the same way, the backtracking with residual agrees with the forward tracking without residual under the condition |wu| + |wv| + P < w. To make both methods compatible, the condition is reduced as |wu| + |wv| + max(E, P ) < w.
According to this statement, what really matters is the evaporation from source A at time t 1 and the precipitation at B at time t 2 . Values such as evaporation in other time or precipitation in other places are not important. As shown in Fig. 4, the precipitation can be any value (under certain conditions described later), as long as-for example-the precipitation in Shandong from 18 to 20 July remains unchanged, the contribution of region 2 stays the same. The highprecipitation center around 32 • N, 110 • E in Fig. 1 actually has no influence on region 2 to its south and the source-receptor relation in Shandong to its east.
The accuracy of the results needs further discussion. First, there is a basic assumption in both models: the assumption of well-mixed atmosphere. In practice, the well-mixed assumption may not be satisfied. Fitzmaurice (2007) pointed out three categories of precipitation that correspond to different mix conditions. One is convection, where the precipitation has a local bias since local moisture contributes more than its proportion. The second is upper troposphere storms, where the advected moisture contributes more. The last is frequent deep convection, where the original well-mixed assumption is valid, such as the rainfalls in Thailand in the monsoonal period (Yoshimura et al., 2004). Several methods have been proposed to account for the non-well-mixed condition (e.g., Burde, 2006;Fitzmaurice, 2007), but they depend on empirical parameters that are seldom available and are not available on transient timescales. Van der Ent et al. (2013) divided their original WAM into a two-layer model to account for this condition. Within each layer, it is well-mixed. The results match well with a regional climate model-based moisture tracking method, which is taken for reference. Second, the unclosed nature of reanalysis data introduces uncertainties. In this study, the GPCP product serves as another precipitation dataset that provides more reliable spatial distribution of rainfall. The resulting contribution map derived from GPCP data is only slightly different than that from the ERA-I data, with the exception that the GPCP map is a little higher and the range of the same contribution level expands more. Besides precipitation, evaporation is a less reliable variable. As analyzed before, the source region's evaporation and the sink region's precipitation are crucial volumes that directly influence the results.
Third, different schemes for dealing with the residual also influence the final result. The backward tracking with residual brings more spatial variations than changing a precipitation dataset in this study, because changing the precipitation only influences the magnitude of the result; the spatial pattern basically stays the same. The residual is more spatially variable, which brings more spatial variations to the result.
Due to many uncertainties, the results we obtained from either forward tracking or backward tracking are more suitable for reference than for accuracy. The spatial pattern of recycling moisture is more reliable than the numerical value.
In this paper, the moisture sources for the selected rainfall event are concentrated to the southwest of Shandong as in Figs. 2 and 3, which are of terrestrial origin. The study of Chen et al. (2011), who investigated moisture sources in a similar case using a Lagrangian method, found that the terrestrial evaporation is more important than oceanic source. Of the terrestrial evaporation, they found that the Indo-China Peninsula, Sichuan, and Yunnan provinces in China contributed more. Their result is consistent with the result in this paper.
For summer precipitation in an extratropical area like Shandong Province of China, the precipitation might be of convection type. In this case, the precipitation has a local bias. As in Fig. 2, the nearer to the precipitation sink upwind, the higher the proportion that the evaporation there becomes involved in the precipitation process; the further away, the lower the proportion. Thus, there is an aggregated effect in the upwind area near the precipitation sink, where the area's moisture contribution counts more than calculated and dispersed effect in regions that are further away. Similarly, there is a more dispersed effect when the tagged water precipitates downwind far away from the source region (as in Fig. 4).

Conclusions
This paper applies a modification WAM to track the sources of an extreme precipitation event in Shandong Province of China from 18 to 20 July 2007. In the experiments of forward tracking and backward tracking, we found that backward tracking without the residual agreed with forward tracking with consideration of the residual. This result is validated by experiments conducted in this case study and other cases not shown here. We found that there is a mathematical necessity in the difference schemes, so we provided key demonstration steps and the conditions under which the two results completely matched. In addition, we also found that forward tracking without the residual agreed with backward tracking with the residual under the same conditions. Backward moisture tracking (without residual) of the Shandong precipitation event in July 2007 (nearly 58.7 mm in 3 days, based on the ERA-I data) indicated that the three source regions contribute (from near to far) 7.8, 4.3, and 1.4 mm of water, which account for 13.3%, 7.3%, and 2.4%, respectively. Using the GPCP data (63.2 mm in total), the source regions contribute 7.8, 4.5, and 1.5 mm, which account for 12.3%, 7.2%, and 2.4%, respectively. The forward tracking method with residual precisely matches these results. In backward moisture tracking (with residual) based on ERA-I data, the results are 6.9, 3.5, and 0.7 mm, which are not the same as backward moisture tracking (without the residual). The results of forward tracking (without the residual) precisely match the results of backward tracking (with the residual).
Due to uncertainties in data, different schemes for dealing with the residual, and the well-mixed assumption, the moisture tracking result is more suitable for reference. However, the spatial pattern is more reli-able. In any case, we are able to determine that the moisture in this Shandong precipitation event originated mostly from the near upwind area of Southwest China, which is of terrestrial origin, and the neighboring West Pacific contributed little.