How can the zero-lift drag coefficient (parasitic drag) be calculated?
Yes, the wetted area is roughly twice the reference area. Now details depend on how well the reference area captures the exposed area of the wing - dihedral will already increase the wetted area by a factor proportional to the inverse of the cosine of the dihedral angle.

But there is more. Airfoil thickness means that the air has to flow around the airfoil. This displacement effect causes the flow around a thick airfoil to speed up more than around an equivalent but thinner airfoil. The thicker airfoil pushes the air aside and around itself more, causing the flow to accelerate and create more friction than the slower flow around a thinner airfoil. This effect is normally approximated with an additional term in the friction drag formula which is proportional to relative thickness.

Next, the type of boundary layer flow needs to be known. Rough surfaces or high sweep angles will provoke an early transition from laminar to turbulent flow. Read this answer for a more detailed discussion.

Another correction is needed for the Mach number, even in subsonic flow. Of course, once flow becomes transsonic or supersonic, wave drag needs to be added, too.

First you need to calculate the friction coefficient which depends on the Reynolds and Mach numbers of your airfoil flow and the relative mean roughness R:
$$c_f = \frac{\frac{0.43}{log(100/R)^{2.56}}-\frac{1700}{100/R}}{\sqrt{1+0.14\cdot Ma^2}}$$

Next, you approximate airfoil drag as explained above:
$$c_{d0} = c_f\cdot \left(2 + 4\cdot\delta + 120\cdot\left(\frac{1}{\sqrt{1-Ma^2}}\right)^3\cdot\delta^4 - 0.09\cdot Ma^2\right)$$
where $\delta$ is the relative thickness of your airfoil.

The term $\frac{1700}{100/R}$ in the friction drag equation allows for the initially laminar boundary layer. Change the factor 1700 depending on how much laminarity your airfoil offers. This answer shows a graph with the possible range. In the airfoil zero lift drag formula you see first the factor 2 which accounts for the fact that the wing has two sides. To that you add the thickness summand to allow for the displacement effect. The third term with the Prandtl-Glauert factor shows that the formula only works well for Mach &lt; 1, and both the third and fourth term are empirical factors to improve accuracy over Mach.