Swapping your use of $k$ and $l$ in the second term gives
$$\sigma(n) = \sum_{k=1}^{n-1}\frac{\sigma((n,k))}{n-1}+\frac{n}{n-1}\sum_{k=1}^{n-1}\sum_{l=1}^{(n,k)}\frac{\sigma((n,k,l))}{(n,k)}$$
Multiplying by $n-1$,
$$(n-1)\sigma(n)=\sum_{k=1}^{n-1}\left(\sigma((n,k))+n\sum_{l=1}^{(n,k)}\frac{\sigma((n,k,l))}{(n,k)}\right)$$
By taking the sum up to $k=n$ rather than $k=n-1$ we get
$$n\sigma(n)+\sum_{l=1}^{n}\sigma((n,l))=\sum_{k=1}^{n}\left(\sigma((n,k))+n\sum_{l=1}^{(n,k)}\frac{\sigma((n,k,l))}{(n,k)}\right)$$
Then after getting rid of the redundant terms and dividing both sides by $n$ we have
$$\sigma(n)=\sum_{k=1}^{n}\sum_{l=1}^{(n,k)}\dfrac{\sigma((n,k,l))}{(n,k)}$$
The function $a(n)$ is multiplicative. If we define
\begin{align*}
f(m) &= \begin{cases}
\sqrt m, &\text{if $m$ is a perfect square}, \\ 0, &\text{otherwise},
\end{cases} \\
g(m) &= \frac{\sigma^*(m)}{2^{\omega(m)}},
\end{align*}
then it's reasonably easy to see that both $f$ and $g$ are multiplicative functions and that
$$
a(n) = \sum_{k\mid n} f(k) g(n/k).
$$
In other words $a=f*g$ is a Dirichlet convolution of two multiplicative functions and is therefore multiplicative. Its value at the prime power $p^j$ is
$$
\sum_{d^2\mid p^j} d \frac{\sigma^*(n/d^2)}{2^{\omega(n/d^2)}} = \sum_{i=0}^{\lfloor j/2\rfloor} p^i \frac{\sigma^*(p^{j-2i})}{2^{\omega(p^{j-2i})}}.
$$
If $j$ is odd then this is
$$
\frac12 \sum_{i=0}^{\lfloor j/2\rfloor} p^i \sigma^*(p^{j-2i}) = \frac12 \sum_{i=0}^{\lfloor j/2\rfloor} p^i (p^{j-2i}+1) = \frac12\sum_{i=0}^{(j-1)/2} (p^{j-i}+p^i),
$$
which can be evaluated via finite geometric series; a similar calculation will cover the case where $j$ is even (careful with the $2^\omega$ term when $i=\frac j2$).
Best Answer
The convolution of two multiplicative functions is multiplicative. Since $\sigma=id\ast 1$, we have $\sigma$ is multiplicative. Even more, we have $f(n)=\sum_{d|n}\sigma(d)$ is multiplicative($f=\sigma\ast 1$, the convolution of two multiplicative functions).
Multiplicative means that $f(mn)=f(m)f(n)$ when $\gcd(m,n)=1$. In particular a multiplicative function is determined by the values of it at prime powers.
In particular, we have $2020=(2^2)(5)(101)$. Hence, $$f(2020)=f(2^2)f(5)f(101)=(\sigma(1)+\sigma(2)+\sigma(2^2))(\sigma(1)+\sigma(5))(\sigma(1)+\sigma(101))=(1+3+7)(1+6)(1+102)=7931$$