All $\mathbb{Z}$-submodules of $\mathbb{Z}\oplus\mathbb{Z}$ are free. 
We are looking for a somewhat direct proof that all $\mathbb{Z}$-submodules of $\mathbb{Z}\oplus\mathbb{Z}$ are free.

We already know that if $A$ is a P.I.D. and we have a free $A$-module $M$, then every sub $A$-module of $M$ is free. However, the proof for this requires quite advanced arguments and we have a strong feeling that a simple proof can be provided for our specific case.
We tried using projections but did not manage to write a coherent proof of the statement.
Thank you very much in advance for your tips and tricks.
 A: Thanks to @Hwang for the idea on this proof. Here it goes.
Suppose $M$ is a sub-$\mathbb{Z}$-module of $\mathbb{Z}^2$. We take $v_1=(p,0)\in M$ such as $p$ is positive and minimal. If such a $p$ does not exist, we can simply not include this element in our basis.
Now choose $v_2=(q,r)\in M$ for $r$ positive and minimal. Again, if such a point does not exist, it means that all points are on the axis and we can take $v_1$ as our basis.
First of all, we have $\langle\{v_1,v_2\}\rangle\subseteq M$, since $M$ is a module and $v_1,v_2\in M$.
Now, choose $(a,b)\in M$. We show that it can be expressed as a linear combination of $v_1$ and $v_2$. If $b\neq 0$, it means that a choice for $r$ exists and thus $r\neq 0$. We must have that $r|b$. Indeed, if $r$ does not divide $b$, we have that $d=\gcd(r,b)<r$ and by linear combination of $(a,b)$ and $(q,r)$ we could by Bézout's identity find an element for which the 2nd coordinate is $d<r$, a contradiction with our choice of $r$.
Since $r|b$, we have $b=kr$. We can then look at $(a,b)-k(q,r)=(a-kq,0)\in M$. Now, by a similar argument as before, if $a-kq\neq 0$ (and thus $p\neq0$), we can find that $p|(a-kq)$, say $pt=a-kq$. This means we have found that $(a,b)=t v_1+k v_2$.
Now, for linear independance, if we had that no such $p$ or $r$ existed, we excluded respectively $v_1$ or $v_2$ from the basis, so we can assume both are in the basis. Then, if we have a linear combination $k_1(p,0)+k_2(q,r)=(0,0)$, we can deduce $k_2=0$, then $k_1=0$.
A: A $\mathbb{Z}$-module is just an abelian group. Since any submodule $M\subset \mathbb{Z}^2$ has rank $\dim_{\mathbb{Q}} M\otimes \mathbb{Q} \leq 2$, it is finitely generated. It's clearly torsion-free, so it's free as a $\mathbb{Z}$-module.
A: How about finding a basis explicitly? Let $M$ be a submodule. Pick primitive $v_1 \in M$. We may assume $v_1 = (p,0)$. Pick $v_2 = (q,r) \in M$ with smallest $r>0$ if such element exists. We can check $\{v_1, v_2\}$ is a basis for $M$. Write an element of $M$ as a linear combination of $v_1$ and $v_2$ and show coefficients are integers.
A: I can suggest a slight simplification, coming from the fact that $\mathbb Z$ is a Euclidean domain (not just a principal ideal domain) - though I would be very pleased if somebody could suggest a more elementary answer than this!
Suppose $M'$ is a submodule of $\mathbb Z^{\oplus m}$. (In your example, $m = 2$.) Then $M'$ is finitely generated as a $\mathbb Z$-module. (For example, $\mathbb Z$ is a PID, hence $\mathbb Z$ is noetherian, hence $\mathbb Z^{\oplus m}$ is noetherian, hence all submodules of $\mathbb Z^{\oplus m}$ are finitely generated...)
So let $x_1, \dots, x_n$ be a finite set of generators for $M'$. We can define a homomorphism $f : \mathbb Z^{\oplus n} \to \mathbb Z^{\oplus m}$ by sending $(a_1, \dots, a_n)\mapsto \sum_i a_i x_i$.
Since $\mathbb Z$ is a Euclidean domain, it is possible to choose bases for $\mathbb Z^{\oplus n}$ and $\mathbb Z^{\oplus m}$ in which $f$ is represented by a diagonal matrix. This is a consequence of the Smith normal form algorithm. Therefore, $M'$ is a free $\mathbb Z$-module, whose rank is equal to the number of non-zero diagonal entries in the matrix.
