## 3 Methodology  
We analyzed the ship flow sequences of three representative container ports by the MF-DFA method, which is a generalization of DFA method [35; 36]. To obtain the generalized Hurst exponent, we followed the five-step procedure introduced by [9; 36] to measure the multifractality and nonstationary behavior of Brazilian rivers. [PERSON] et al. [9] pointed out that the periodic components in the sequence should be removed at the first stage before beginning with the general procedure of [PERSON] et al. [36].  
For a record \(x(i),i=1,2,...,N\), where \(N\) denotes the length of the record, the MF-DFA consists of the following steps [35; 36]:  
Step 1: We first integrate the series and obtain the profile \(y(j)\),  
\[y(j)=\sum_{i=1}^{j}[x_{i}-\mu], \tag{1}\]  
where \(\mu\) is the mean value of the entire series.  
Step 2: The integrated series \(y(j)\) is divided into boxes of equal length \(s\).  
Step 3: In each box of length \(s\), we calculate a polynomial fitting of \(y(j)\), which represents the trend in that box. The shape of the polynomial trend is defined by the order \(m\). A higher order \(m\) yield a more complex shape of the trend, but might lead to overfitting for a time series within small segment sizes. Therefore, in this study, we choose \(m\) equal to 2 as suggested by [PERSON] (2012) [37]. The \(y\) coordinate of the fit line in each box is denoted by \(y_{s}(j)\).  
Step 4: The integrated series \(y(j)\) is detrended by subtracting the local trend \(y_{s}(j)\) in each box of length \(s\).  
Step 5: For a given box size,  
\[F_{q}(s)=\Bigg{[}\frac{1}{N}\sum_{j=1}^{N}|y(j)-y_{s}(j)|^{q}\Bigg{]}^{\frac{ 1}{q}}. \tag{2}\]  
If the series is long-range and the power-law is correlated for large values of the time scales, the fluctuation functions \(F_{q}(s)\) can be written as Equation (3):  
\[F_{q}(s)\sim s^{H(q)}, \tag{3}\]  
where \(H(q)\) is the generalized Hurst exponent. The generalized Hurst exponent \(H(q)\) can be obtained by observing the slope of the log-log plots of \(F_{q}(s)\) and scale \(s\) through the method of ordinary least squares (OLS). If \(H(q)\) is independent of \(q\), the time series is monofractal, otherwise, it is multifractal. Table 1 shows the relationship between the long-range dependence of time series and the generalized Hurst exponent.  
The value of the generalized Hurst exponent equal to 0.5 indicates an uncorrelated time series. A generalized Hurst exponent value larger than 0.5 indicates a positive long-range dependence and persistence of the series. In other words, the larger the \(H\) value is, the stronger the persistence is. A generalized Hurst exponent value smaller than 0.5 indicates a negative long-range dependence and antipersistence of the series. This means the closer the \(H\) value is to 0, the stronger the antipersistence is.  
When \(q\) is equal to two, \(H(2)\) is identical to the well-known Hurst exponent. Generally, the Hurst exponent is between 0 and 1. However, it is worth noting that the generalized Hurst index obtained by applying the MF-DFA method in this study may be greater than 1 [36; 38].  
The singularity spectrum \(f(\alpha)\) is introduced to measure the degree of multifractality of the series and can be obtained through Legendre transform:  
\[\alpha=H(q)+qH^{\prime}(q), \tag{4}\]  
\[f(\alpha)=q[\alpha-H(q)]+1. \tag{5}\]  
Here, \(\alpha\) is the singularity strength and used to characterize the singularities of the time series. \(f(\alpha)\) indicates the dimension of the subset of sequences that is characterized by \(\alpha\). The strength of multifractality can be estimated by the spans of singularity given by  
\[\Delta\alpha=\alpha_{\text{max}}-\alpha_{\text{min}}. \tag{6}\]