## 2 Methods
### Experimental Designs  
We designed the experiments based on the RF tuning parameters, whose values are listed in Table 1. Each parameter is defined as follows: \(L\), the horizontal correlation length scale, determines how much perturbed errors propagate in a horizontal direction; \(\tau\), the de-correlation time scale, determines how long the perturbed errors will be sustained; \(\sigma\), the standard deviation, controls the amplitudes of RF; and \(\mu\), the tapering function, determines whether the error would exponentially decrease or remain the same in the vertical layers. With these parameter settings, each stochastic perturbation is activated simultaneously in a stable state. Under the same RF, the tendency fluctuation in SPDT is larger than that in SPPT; thus, in order to suppress spurious instability, we designed the SPDT to generate sufficiently small RF by setting smaller \(\sigma\) while activating \(\mu\) to render the RF smaller with height [25].  
We conducted two sets of EDA runs--the control run (CTL) and the stochastic perturbation run (STO)--to assess the SPHT scheme. The additive inflation and the RTPS were included in both runs; however, the SPDT and the SPPT were simultaneously activated only in the STO. Thus, the differences between the results from CTL and STO are solely due to the stochastic perturbation. The experiments started at 12:00 UTC on 22 June 2018 and ended at 18:00 UTC on 7 July 2018; this period was long enough to examine the response of BEC in an EDA system. We specified the first 78 h as a spin-up period and excluded this period from the error analyses.  
We assessed the proposed scheme by quantifying the ensemble features. The ES was examined with regard to the areas that were underestimated or overestimated, compared to the EME that is defined as the root-mean-square deviation (RMSD) with respect to the ensemble mean. The EME was evaluated against the Integrated Forecast System (IFS) analysis from the European Centre for Medium-Range Weather Forecast (ECMWF), having a horizontal resolution of 25 km with 25 pressure levels from 1000 to 1 hPa. Currently, the analysis of IFS/ECMWF shows the smallest RMSD compared with the analyses of other operational centers (see [42]); thus, we regard the IFS analysis as the true state against which the EME is estimated.We also checked the performance of EDA by the new SPHT scheme through the _accuracy_ and _reliability_. The accuracy was measured by the reduction of EME  
\begin{table}
\begin{tabular}{c c c c c} \hline \hline
**Experiment** & \(L\) & \(\tau\) & \(\sigma\) & \(\mu\) \\ \hline
**SPDT** & 500 & 10,800 & 0.5 & on \\
**SPPT** & 500 & 21,600 & 1.0 & off \\ \hline \hline \end{tabular}
\end{table}
Table 1: Configuration of the random forcing (RF) tuning parameters for the experiment using stochastically perturbed dynamical tendency (SPDT) and stochastically perturbed physical tendency (SPPT): \(L\) is the horizontal correlation length scale (in km), \(\tau\) the de-correlation time scale (in s), \(\sigma\) the standard deviation, and \(\mu\) the tapering function.  
while the reliability was evaluated by comparing the ES and the EME [43]. In particular, we define the reliability index (\(\gamma\)) as  
\[\gamma=\frac{\text{ES}}{\text{EME}}, \tag{3}\]  
which depicts the balance between the ES and the EME.