Step 1: Formula for surface area.
For a surface \( z = f(x, y) \), the area is
\[
A = \iint_D \sqrt{1 + \left( \frac{\partial z}{\partial x} \right)^2 + \left( \frac{\partial z}{\partial y} \right)^2} \, dx \, dy.
\]
Given \( z = 16 - x^2 - y^2 \),
\[
\frac{\partial z}{\partial x} = -2x, \frac{\partial z}{\partial y} = -2y.
\]
So,
\[
A = \iint_D \sqrt{1 + 4(x^2 + y^2)} \, dx \, dy.
\]
Step 2: Convert to polar coordinates.
\[
x = r \cos \theta, \; y = r \sin \theta, \; dx \, dy = r \, dr \, d\theta.
\]
Since \( z = 0 \) corresponds to \( r^2 = 16 \),
\[
A = \int_0^{2\pi} \int_0^4 \sqrt{1 + 4r^2} \, r \, dr \, d\theta.
\]
Final Answer: \[ \boxed{A = \int_0^{2\pi} \int_0^4 \sqrt{1 + 4r^2} \, r \, dr \, d\theta.} \]