## 1 Introduction  
Numerical weather prediction (NWP) includes inevitable forecast errors due to uncertainties in both initial conditions and models. The model error would be observed even though the initial conditions perfectly describe the true state of atmospheric flows [1]; thus, the differences in tendencies between the model and the real world always exist. As an alternative to deterministic forecasting, which has limitations in representing this uncertainty, an ensemble prediction system (EPS) can be utilized to overcome the limitation [2, 3]. In particular, the ensemble data assimilation (EDA) system is beneficial to represent the initial uncertainties and flow-dependent background error covariance (BEC), thus improving the assimilation of observations [4]. If the ensemble perturbation includes the true observation and the model error, the unperturbed analysis will represent the analysis error [5]. However, the EDA usually suffers from the underestimation of BEC due to the limited ensemble size, sampling errors, and model error [6]. This underestimation leads to a filter divergence problem: the analysis diverges from the natural state by ignoring the observation influence due to a small forecast uncertainty. Conversely, with an over-dispersive ensemble, the analysis ignores the model errors due to a large forecast uncertainty. If plenty of ensemble sizes are available, the sampling error is expected to be moderately eliminated [7,8]; however, this would be impractical due to limited computational resources. In general, the problems related to covariance underestimation and sampling error can be solved with covariance inflation and localization, respectively [6]--the former inflates the small uncertainty while the latter reduces sampling errors due to limited ensemble size. In this study, we focus on the covariance inflation only. During the last a few decades, several inflation methods have been proposed, including multiplicative inflation (e.g., [9]), additive inflation (e.g., [10,11]), and some alternative inflation methods using relaxation (e.g., [12, 13, 14]) and adaptive approaches (e.g., [6, 15, 16, 17]). Nevertheless, it has been found that the ensemble BEC still remains under-dispersive and underestimates the true uncertainty [18, 19]. In addition to imperfect initial conditions, forecast errors are attributed to model uncertainties, stemming from approximations in governing equations, truncation and round-off errors, misrepresentations of physical processes, etc. The stochastic approach with random forcing (RF) is widely used to represent the model uncertainties in various applications [20] as well as in the EPSs of many operational centers [1]. The stochastically perturbed physical tendency (SPPT) scheme has been proven to be effective in ensemble spread (ES) and ensemble mean error (EME), compared to the traditional inflation methods [4, 21]. Here, the ES--a standard deviation from the ensemble mean--provides an estimate of the forecast uncertainty [22], and its average is expected to have a similar spread to the EME. Recently, several studies have focused on the dynamic part of a forecast model; for instance, the stochastic kinetic energy backscatter (SKEB) assumes that the model error is associated with inadequate scale interactions in truncated NWP models [1, 23, 24]. Although it showed remarkable performance in the medium-range EPS, it has been difficult to obtain a positive impact in EDA [11]. As an alternative method, [PERSON] and [PERSON] [25] proposed the stochastically perturbed dynamical tendency (SPDT) scheme by assuming that the model error is associated with computational representations of the underlying partial differential equations to solve the atmospheric motion. They reported a significant improvement in prediction skills in a seasonal ensemble framework by perturbing both dynamical and physical tendencies, especially in poleward transient eddy and seasonal precipitation, compared to perturbing the physical tendency only.