For thermal dissipation sap flow time series

The variables listed in the table below are mandatory for each method to determine zero-flow conditions (SP: Successive predawn method; PD: Daily predawn method; MW: Moving window method; DR: Double regression method; ED: Environmental dependent method).

\(\Delta T\) is the temperature difference between sap flow probes. You can select the multiple methods at once during the automatic processing.

Here are the more details of each variable column.

• Timestamp: A timestamp of class POSIXct or POSIXt, indicating the timings of the end of each measurement in local time. Any interval (typically 15 to 60 min) is allowed, but the timestamps must be equally spaced and arranged chronologically.

• \(\Delta T\): A numeric value of the temperature difference between sap flow probes. The units should be in \(^\circ \mathrm{C}\). It is acceptable to include missing values, but the label representing the missing value must be unique such as -9999.

• Global solar radiation: A numeric value indicating incident short-wave radiation. The units should be in W m^-2. If you choose the PD, MW, DR, or ED methods, the time series must be included and gap-filled previously. These methods use it only to define daytime, so you can use PPFD (photosynthetic photon flux density) or another variable with a strong correlation with global solar radiation. In such cases, you must modify the threshold value for daytime definition specified by the argument thres_radi in run_fluxfixer() or calc_dtmax(). In case you have no radiative time series, this package offers the calc_sw_in_toa() function that estimates the global solar radiation at the top of the atmosphere at a specific location. You can add this time series to your input data frame and determine the zero-flow conditions. See more details in ?calc_sw_in_toa().

• Air temperature: A numeric value and its units should be in\(^\circ \mathrm{C}\). If you select the ED method, this time series must be included and gap-filled previously.

• Vapor pressure deficit (VPD): A numeric value representing the difference between saturated water vapor pressure and the actual partial water vapor pressure in the atmosphere. The units should be in hPa. If you select the ED method, this time series must be included and gap-filled previously. And even if you do not select the ED method, the VPD is recommended to be included in the input data frame to construct a more powerful random forest model for outlier detection and gap-filling processes, as the VPD typically controls the opening of plant stomata. In this case, you do not have to fill gaps in the time series previously.

• Other variables: Even when you select only the SP method, you need to add as many drivers that influence sap flow dynamics as possible in the input data frame to construct a more powerful random forest model for outlier detection and gap-filling processes. For example, if you want to post-process sap flow data observed in a deciduous forest, you are recommended to include a time series representing forest leaf amount, such as LAI (leaf area index).

You can also include columns that are unrelated to the post-processing, but be sure to indicate which variables are to be used in the random forest construction when you execute run_fluxfixer() by specifying vctr_colname_feature.

For time series unrelated to thermal dissipation sap flow

In this case, the input data frame is required to have the columns below.

• Timestamp: A timestamp of class POSIXct or POSIXt, indicating the timings of the end of each measurement in local time. Any interval is allowed, but the timestamps must be equally spaced and arranged chronologically.

• Target variable: A numeric value of the variable you want to post-process. It is acceptable to include missing values, but the label representing the missing value must be unique such as -9999.

• Other variables: You need to add as many drivers that influence targeted variable dynamics as possible in the input data frame to construct a more powerful random forest model for outlier detection and gap-filling processes.

You can also include columns that are unrelated to the post-processing, but be sure to indicate which variables are to be used in the random forest construction when you execute run_fluxfixer() by specifying vctr_colname_feature.

Data Input

For those who jumped into this section directly, please load the sample data first.

## Load sample data
data("dt_noisy")

Let’s explore the functions deeply that comprise the automated function, run_fluxfixer().

Anomaly correction

To begin with, you remove anomalies and correct the time series. This process is divided into three procedures: the absolute limit test, the short-term drift correction, and the high-frequency noise removal.

## Conduct the absolute limit test
df_all <-
  dt_noisy %>% 
  dplyr::mutate(dt_al = check_absolute_limits(vctr_target = dt))

## Specify the period of the short-term drift
time_drft_head <- as.POSIXct("2013/05/14 13:30", tz = "Etc/GMT-8")
time_drft_tail <- as.POSIXct("2013/05/17 15:00", tz = "Etc/GMT-8")

## Conduct the short-term drift correction
df_all <-
  df_all %>% 
  dplyr::mutate(dt_sdc = modify_short_drift(vctr_time = time,
                                            vctr_target = dt_al,
                                            vctr_time_drft_head = time_drft_head,
                                            vctr_time_drft_tail = time_drft_tail))

Outlier removal

The next step is to remove outliers in the target time series. This process consists of two procedures: Z-score outlier removal and random forest outlier removal.

## Specify the sensor replacement timing
time_prd_tail <- as.POSIXct("2013/05/14 13:00", tz = "Etc/GMT-8")

## Conduct the Z-score outlier removal
df_z <-
  remove_zscore_outlier(vctr_time = df_all$time,
                        vctr_target = df_all$dt_sdc,
                        vctr_time_prd_tail = time_prd_tail)

df_all <-
  df_z %>% 
  dplyr::select(z_cleaned, avg_cleaned, sd_cleaned) %>% 
  dplyr::bind_cols(df_all, .)

## Conduct the random forest outlier removal
df_rf <-
  remove_rf_outlier(df = df_all,
                    colname_label = "z_cleaned",
                    vctr_colname_feature = c("p", "sw_in", "ta", "vpd", "ws", "swc"))

df_all <-
  df_all %>% 
  dplyr::mutate(z_rf = df_rf$stats$cleaned)

## Calculate the reference values of average and standard deviation
vctr_stats_ref <-
        data.frame(z_target = df_all$z_rf,
                   avg_target = df_z$avg_cleaned,
                   sd_target = df_z$sd_cleaned) %>%
        dplyr::mutate(target = ifelse(z_target != -9999 &
                                        avg_target != -9999 &
                                        sd_target != -9999,
                                      z_target * sd_target + avg_target,
                                      -9999)) %>%
        dplyr::pull(target) %>%
        calc_ref_stats(vctr_time = df_all$time,
                       vctr_target = .,
                       vctr_time_prd_tail = time_prd_tail)

Gap-filling

fill_gaps() replaces all missing values in a target time series with values estimated by a random forest model. The hyperparameter calibration process is the same as that of the random forest outlier removal. After determining the optimal hyperparameters, they are used to construct the optimal random forest model. Predicted time series are equal to average from 500 (default; specified by n_tree) tree outputs at each time point. If the input targeted value is missing, the predicted value is used for the imputation.

## Fill all the missing values by the random forest model
df_gf <-
  fill_gaps(df = df_all,
            colname_label = "z_rf",
            vctr_colname_feature = c("p", "sw_in", "ta", "vpd", "ws", "swc"))

df_all <-
  df_all %>% 
  dplyr::mutate(z_gf = df_gf$stats$gapfilled)

The column z_gf in the df_all data frame represents the gap-filled standardized \(\Delta T\) time series.

Time series retrieval

retrieve_ts() converts the standardized Z-score time series into a time series in its original units. This retrieval is performed by multiplying the Z-score by the standard deviation and adding the average. If the average and standard deviation time series are the same as those in converting the original time series into the Z-score time series, the original values with the original average and standard deviation are retrieved. If reference values of the average and/or standard deviation are used, the output time series are detrended and/or applied to signal damping correction.

For the sample data process, detrending is needed. So, you set detrend as TRUE, input the reference average calculated by calc_ref_stats() into avg_ref, and let correct_damping remain FALSE.

## Retrieve the time series in its original units
df_all <-
  df_all %>%
  dplyr::mutate(dt_processed = retrieve_ts(vctr_target_z = df_all$z_gf,
                                           vctr_target_sd = df_all$sd_cleaned,
                                           detrend = TRUE,
                                           avg_ref = vctr_stats_ref[1]))

The column dt_processed represents the quality-controlled \(\Delta T\) time series, and the same as the output time series by run_fluxfixer() in the previous section. Well done!

In the following parts, you will learn the processes specific to thermal dissipation sap flow data.

Zero-flow conditions determination

calc_dtmax() provides multiple \(\Delta T_\mathrm{max}\) time series estimated by different methods below:

• The successive predawn (SP) method defines the \(\Delta T_\mathrm{max}\) for a day as the maximum \(\Delta T\) within a 24-hour period that begins at 5:00 a.m. (default; just before daybreak in temperate zones and specified by thres_hour_sp). In other words, the day starts at predawn, not midnight, and the maximum value for that period is assumed to be \(\Delta T_\mathrm{max}\). This method has the advantage of being able to calculate \(\Delta T_\mathrm{max}\) quickly while minimizing the effect of nocturnal transpiration on \(\Delta T_\mathrm{max}\) estimation.

• The daily predawn (PD) method defines the \(\Delta T_\mathrm{max}\) for a day as the maximum \(\Delta T\) between midnight and the morning (8:00 a.m. in local time, specified by thres_hour_pd) when the global solar radiation is below the threshold value, specified by thres_radi. See more details in Peters et al. (2018).

• The moving window (MW) method selects the maximum value of \(\Delta T_\mathrm{max}\), estimated by the PD method, using a moving window with an eleven-day length (default; specified by wndw_size_dtmax). The selected \(\Delta T_\mathrm{max}\) is considered to be the final \(\Delta T_\mathrm{max}\). See more details in Peters et al. (2018).

• The double regression (DR) method first calculates the moving window mean value of \(\Delta T_\mathrm{max}\), estimated by the PD method, with an eleven-day length (default; specified by wndw_size_dtmax). The \(\Delta T_\mathrm{max}\) that is lower than the mean is omitted, and then the moving window mean is recalculated as the final \(\Delta T_\mathrm{max}\). See more details in Peters et al. (2018).

• The environmental dependent (ED) method filters the \(\Delta T_\mathrm{max}\), estimated by the PD method, using the environmental conditions when plants let their sap flow nearly zero. A stable \(\Delta T\), with a low coefficient of variation, and low air temperature or vapor pressure deficit over a two-hour period, characterizes these zero-flow conditions. The threshold values can be changed by thres_cv, thres_ta, and thres_vpd arguments. See more details in Oishi et al. (2016) and Peters et al. (2018). After the filtering, the daily \(\Delta T_\mathrm{max}\) is interpolated if necessary.

By default, calc_dtmax() only outputs \(\Delta T_\mathrm{max}\) time series by the SP method. You can select any combination of the methods by specifying method, such as method = c("sp", "pd"). Make sure you indicate the vectors of the meteorological variables required for the selected methods into vctr_radi, vctr_ta, and vctr_vpd arguments. See the table in the “Data preparation” section to check which variables are needed for each method.

## Calculate dTmax by the SP and PD methods
df_dtmax <-
  calc_dtmax(vctr_time = df_all$time,
             vctr_dt = df_all$dt_processed,
             vctr_radi = df_all$sw_in,
             method = c("sp", "pd"))

df_all <-
  df_dtmax %>% 
  dplyr::select(dtmax_sp, dtmax_pd) %>% 
  dplyr::bind_cols(df_all, .)

The columns dtmax_sp and dtmax_pd represent the \(\Delta T_\mathrm{max}\) time series calculated by the SP and PD method, respectively.

Sap flux density calculation

This is the final step. calc_fd() calculates \(F_\mathrm{d}\) (sap flux density) time series, including heartwood correction. \(F_\mathrm{d}\) is estimated using a power-type function introduced by Granier (1985, 1987). First, a dimensionless index K is obtained from \(\Delta T\) and \(\Delta T_\mathrm{max}\). Second, K is raised to the power beta and then multiplied by alpha, obtaining \(F_\mathrm{d}\).

If the sapwood width is shorter than the probe insertion length, \(\Delta T\) can be overestimated, resulting in an underestimation of \(F_\mathrm{d}\). Therefore, heartwood correction is required to correct \(\Delta T\). Optionally, before calculating \(F_\mathrm{d}\), \(\Delta T\) can be replaced with the corrected \(\Delta T\) by specifying the ratio of the probe length to sapwood width. This correction assumes that the \(\Delta T\) measured by the part of the probe that is inserted into the heartwood is always \(\Delta T_\mathrm{max}\). See more details in Clearwater et al. (1999).

For the sample data process, there is no need to do the heartwood correction. So, you let do_heartwood_correction remain FALSE.

## Calculate Fd using the dTmax estimated by the SP and PD methods
df_all <-
  df_all %>% 
  dplyr::mutate(fd_sp = calc_fd(vctr_dt = dt_processed,
                                vctr_dtmax = dtmax_sp),
                fd_pd = calc_fd(vctr_dt = dt_processed,
                                vctr_dtmax = dtmax_pd))

The columns fd_sp and fd_pd represent the calculated \(F_\mathrm{d}\) time series using the \(\Delta T_\mathrm{max}\) estimated by the SP and PD method, respectively.

Congratulations! Now you can use these functions and run_fluxfixer() more appropriately and flexibly. I hope the fluxfixer will make your data quality control tasks less painful.

If you have any questions, please check each function’s help page first, then feel free to contact me.