The Idea of Convolution
My favorite exposition of the topic is in one of Brad Osgood's lectures on the Fourier Transform. The discussion of convolution begins around 36:00, but the whole lecture has additional context that's worth watching.
The basic idea is that, when you define something like the Fourier Transform, rather than working directly with the definition all the time, it's useful to derive higher-level properties that simplify calculations. For example, one such property is that the transform of the sum of two functions is equal to the sum of the transforms, i.e.
$$
F\{f + g\} = F\{f\} + F\{g\}.
$$
That means if you have a function with an unknown transform, and it can be decomposed as a sum of functions with known transforms, you basically get the answer for free.
Now, since we have an identity for the sum of two transforms, it's a natural question to ask what the identity is for the product of two transforms, i.e.
$$
F\{f\}F\{g\} = \space ?.
$$
It turns out that when you calculate the answer, convolution is what appears. The entire derivation is given in the video, and since your question is mostly conceptual, I won't recapitulate it here.
The implication of approaching convolution in this way is that it's an intrinsic part of the way the Laplace Transform (of which the Fourier Tranform is a special case) turns linear constant-coefficient ordinary differential equations (LCCODE's) into algebraic equations. That fact that such a transform is available to make LCCODE's analytically tractable is a large part of the reason why they're studied in signal processing. For example, to quote Oppenheim and Schafer:
Because they are relatively easy to characterize mathematically and
because they can be designed to perform useful signal processing
functions, the class of linear shift-invariant systems will be studied
extensively.
So one answer to the question, is that if you're using transform methods to analyze and/or synthesize LTI systems, sooner or later, convolution will arise (either implicitly or explicitly). Note that this approach to introducing convolution is very standard in the context of differential equations. For example, see this MIT lecture by Arthur Mattuck. Most presentations either present the convolution integral without comment, then derive its properties (i.e. pull it out of a hat), or hem and haw about the strange form of the integral, talk about flipping and dragging, time-reversal, etc., etc.
The reason I like Prof. Osgood's approach is that it avoids all that tsouris, as well as providing, in my opinion, deep insight into how mathematicians probably arrived at the idea in the first place. And I quote:
I said, "Is there a way of combining F and G in the time domain, so
that in the frequency domain the spectrums multiply, the Fourier
transforms multiply?" And the answer is, yes, there is, by this
complicated integral. It's not so obvious. You wouldn't get out of bed
in the morning and write this down, and expect that this was going to
solve that problem. How do we get it? You said, suppose the problem is
solved, see what has to happen, and then we'd have to recognize when
it's time to declare victory. And it's time to declare victory.
Now, being an obnoxious mathematician, you cover your tracks and say,
"Well, I'm simply going to define the convolution of two functions by
this formula."
LTI Systems
In most DSP texts, convolution is usually introduced in a different way (that avoids any reference to transform methods). By expressing an arbitrary input signal $x(n)$ as a sum of scaled and shifted unit impulses,
$$
x(n) = \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}, \tag{1}
$$
where
$$
\delta(n) = \begin{cases}
0, & n \ne 0 \\
1, & n = 0,
\end{cases} \tag{2}
$$
the defining properties of linear time-invariant systems lead directly a convolution sum involving the impulse response $h(n) = L[\space \delta(n) \space]$. If the system defined by an LTI operator $L$ is expressed as $y(n) = L[\space x(n)\space]$, then by applying the repective properties, namely linearity
$$
\overbrace{L[\space ax_1(n) + bx_2(n)\space ]}^\text{Transform of the sum of scaled inputs} = \overbrace{aL[\space x_1(n)\space ] + bL[\space x_2(n)\space ]}^\text{Sum of scaled transforms} , \tag{3}
$$
and time/shift invariance
$$
L[\space x(n) \space] = y(n) \space \xrightarrow{\text{implies}} L[\space x(n-k) \space] = y(n-k), \tag{4}
$$
the system can be rewritten as
$$
y(n) = \overbrace{L\left[ \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}\right]}^\text{Tranform of the sum of scaled inputs} = \overbrace{\sum_{k=-\infty}^\infty{x(k)L[\delta(n - k)]}}^\text{Sum of scaled transforms} = \overbrace{\sum_{k=-\infty}^\infty{x(k)h(n-k) }. }^\text{Convolution with the impulse response}
$$
That's a very standard way to present convolution, and it's a perfectly elegant and useful way to go about it. Similar derivations can be found in Oppenheim and Schafer, Proakis and Manolakis, Rabiner and Gold, and I'm sure many others. Some deeper insight [that goes further than the standard introductions] is given by Dilip in his excellent answer here.
Note, however, that this derivation is somewhat of a magic trick. Taking another look at how the signal is decomposed in $\left(1\right)$, we can see that it's already in the form of a convolution. If
$$
\overbrace{(f * g)(n)}^\text{f convolved with g} = \sum_{k=-\infty}^\infty{f(k)g(n-k)},
$$
then $\left(1\right)$ is just $x * \delta$. Because the delta function is the identity element for convolution, saying any signal can be expressed in that form is a lot like saying any number $n$ can be expressed as $n + 0$ or $n \times 1$. Now, choosing to describe signals that way is brilliant because it leads directly to the idea of an impulse response--it's just that the idea of convolution is already "baked in" to the decomposition of the signal.
From this perspective, convolution is intrinsically related to the idea of a delta function (i.e. it's a binary operation that has the delta function as its identity element). Even without considering its relation to convolution, the description of the signal depends crucially on the idea of the delta function. So the question then becomes, where did we get the idea for the delta function in the first place? As far as I can tell, it goes at least as far back as Fourier's paper on the Analytical Theory of Heat, where it appears implicitly. One source for further information is this paper on Origin and History of Convolution by Alejandro Domínguez.
Now, those are the two of the main approaches to the idea in the context of linear systems theory. One favors analytical insight, and the other favors numerical solution. I think both are useful for a full picture of the importance of convolution. However, in the discrete case, neglecting linear systems entirely, there's a sense in which convolution is a much older idea.
Polynomial Multiplication
One good presentation of the idea that discrete convolution is just polynomial multiplication is given by Gilbert Strang in this lecture starting around 5:46. From that perspective, the idea goes all the way back to the introduction of positional number systems (which represent numbers implicitly as polynomials). Because the Z-transform represents signals as polynomials in z, convolution will arise in that context as well--even if the Z-transform is defined formally as a delay operator without recourse to complex analysis and/or as a special case of the Laplace Transform.